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1.  Introduction 


The  era  of  extracting  gravity  field  information  about  the  earth  from 
only  time-con  Burning  laborious  surface  measurements  is  definitely  over. 

Modern  technology  has  created  highly  sophisticated  instruments  which 
enable  us  to  obtain  different  kinds  of  geodetic  data  with  increasing 
accuracy  at  a speed  never  known  before.  Probably  the  only  technique  cap- 
able of  processing  all  these  data  in  a consistent  way  is  the  method  of 
least-squares  collocation,  which  is  as  beautiful  as  deep  (Morits,  1972). 

The  central  role  in  collocation  is  played  by  the  model  covariance 
function  which  resembles  somehow  the  main  features  of  the  earth's 
gravitational  field.  The  property  of  the  kernel  function  forces  it  to  be 
harmonic  outside  some  internal  sphere;  consequently,  its  eigen-functions 
are  the  spherical  harmonics.  Homogeneity  and  isotropy  postulates 
make  it  dependent  on  only  two  essential  variables:  the  spherical  distance 
between  two  points  and  the  product  of  the  corresponding  geocentric  radii  of 
these  points. 

All  geodetic  data  are  functionals  of  the  disturbing  potential.  Therefore, 
all  covariances  between  geodetically  relevant  quantities  can  be  derived  by 
means  of  the  covariance  propagation  law.  This,  again,  is  the  reason  why 
closed  analytical  expressions  of  the  covariance  function  are  an  absolute 
necessity.  Tscheming  and  Rapp  (1974)  have  derived  closed  expressions 
on  the  basis  of  different  degree  variance  models  and  give  covariance 
expressions  for  geoidal  height,  gravity  anomaly  and  deflection  of  the 
vertical,  most  commonly  used  in  geodetic  applications.  Tscherning  (1976) 
extended  this  work  and  derived  covariance  expressions  for  second 
order  derivatives  of  the  anomalous  potential.  His  very  elegant  and 
convenient  subroutine  COVAX  has  been  widely  applied  by  the  geodetic 
community. 

However,  the  number  of  problems  for  which  the  application  of  COVAX 
becomes  expensive  increases  steadily.  All  these  expensive  problems 
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involve  the  numerical  integration  of  the  covariance  function  resulting  in 
a large  number  of  covariance  computation*.  Two  typical  example*: 
a)  the  prediction  of  mean  gravity  anomalies  from  different  kind*  of 
data  (Rapp,  1977,  1978),  b)  *a tellite- to- staellite  tracking  espec- 

ially in  the  low- low  mode)  (Kryn'ski,  1978).  In  both  case* 
multifold  integrations  of  the  covariance  function  are  involved  (in  the 
first  case  integratione  over  "rectangular"  areas  on  the  sphere,  in  the 
second  case  integrations  along  the  flight  paths  of  two  satellites). 

For  this  reason  the  possibility  had  been  studied  to  cut  down  the 
considerably  the  covariance  computation  time  by  using  approximations  of  the 
covariance  function  (Sunkel,  1978  a,b).  Theoretical  investigations 
indicated  that  this  possibility  exists)  preliminary  practical  calculations 
encouraged  further  detailed  studies.  After  all,  a way  1ms  been  found 
through  the  jungle  which  permitted  a consistent  approximation  of  all 
covariance  expressions  up  to  and  including  second-order  derivatives  of 
the  disturbing  potential.  Basically,  the  method  consists  of  choosing  a 
regular  rectangular  grid  with  respect  to  the  cosine  of  the  spherical 
distance  tf  * co*tj  with  respect  to  the  squared  ratio  s between  the 
radius  of  the  Bjerhammar  sphere  R,  and  some  radius  r > , s,  =(l^  /r,)1, 

and  storing  all  covariances  corresponding  to  these  grid  points  on  some 
file.  All  desired  covariance  expressions  can  be  derived  by  a simple 
differentiation-interpolation  procedure.  The  accuracy  of  the  eo  obtained 
covariancee  can  be  made  arbitrary  high  by  selecting  an  appropriate  small 
grid  pacing. 

This  procedure  provides  the  user  with  the  same  covariances  as 
conveniently  asCOVAX,  but  about  10  times  faster.  The  price  we  pay  is 
mass  storage.  However,  computer  development  statistics  show  clearly 
that  the  amount  of  available  mass  storage  increases  much  more  than  the 
speed-up  of  calculation  time.  Moreover,  since  large  scale  applications 
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of  collocation  (for  which  the  procedure  described  here  ie  designed)  are 
by  their  own  nature  restricted  to  large  computer  systems  with  a large 
central  core  capacity  available,  the  need  for  storage  can  no  longer  be  con- 
sidered an  essential  drawback.  It's  believed  that  this  approximation 
procedure  will  be  especially  valuable  for  collocation  problems  which 
involve  the  integration  of  the  covariance  function. 


2.  The  Spatial  Covariance  Function  and  its  Approximation  Model 

Although  the  choice  of  a particular  covariance  function  model  is  of 
no  concern  to  the  approximation  procedure  (as  long  as  it  is  isotropic),  it 
is  nevertheless  necessary  to  give  a short  description  of  its  essential 
properties  in  order  to  have  a solid  starting  basis. 

The  general  form  of  a homogeneous  and  isotropic  covariance  function 
can  be  expressed  by 

K(P.Q)  = L k,  (iL)n+1  P,  (com,,).  (2-1) 

n=N„  ' / 


where 


P»  Q* . . 


points  in  space, 
geocentric  radii  of  P and  Q, 
spherical  distance  between  P and  Q, 
positive  coefficients, 

Legendre  polynomial  of  degree  n, 

N0. . . starting  value  of  die  summation  (N0>  2), 

R,  . . . radius  of  the  Bjerhammar  sphere. 

K(P,Q)  is  symmetric  with  respect  to  P and  Q,  expressed  by  the  product 
r r*  and  the  dependence  on  cost,q, 


r,  r . . , 

k, ... 
Pt(cosf) 


K(P,Q)  = K(Q,  P); 

it  i a harmonic  outside  the  Bjerhammar  sphere  r = R, , 

A,  K(P.Q)  = Aq  K(P,Q)  = 0. 

An  important  property  which  allows  a two-dimensional  approximation  is 
its  dependence  on  essentially  two  variables  only,  the  spherical  distance 
f,q  and  the  product  r r/.  Since  cost  can  vary  only  between  -1  and  4 1, 
and  Rq/rr'  has  a minimum  value  of  0 for  r -»  • and  a maximum  of  1 
for  r * r#  * R, , the  covariance  functions  domain  is  the  rectangle 


[-1  < t * 1,  0 < s < 1] 


with  the  variables 


t : = cost 


Re 

. • : = -r*/  . 


(2 -2a) 


(2 -2b) 


Operations  on  or  close  to  the  surface  of  the  earth  together  with  locally 
to  regionally  restricted  use,  however,  reduce  the  domain  of  practical 
applications  considerably.  E.  g. : Working  within  a spherical  distance 
range  of  0 s t < 10°  and  within  an  altitude  range  from  0 (earth's  surface) 
to  300  km  (satellite  altitude),  the  operational  domain  reduces  to 

tO.  985  * t * 1,0.999  « s s 0.912J. 

The  following  figure  may  illustrate  the  mapping  given  by  (2-2b): 


Y y . 

y Yy 

r/// 


V////Y 


Fig.  1 (s,  t)  - mapping 


Once  the  operational  domain  hae  been  fixed,  a regular  rectangular 
grid  in  • and  t can  be  arranged.  Two  important  facta  have  to  be  taken 
into  account:  The  numerical  values  of  the  covariances  change  more 
rapidly  for  small  f t he. : for  larger  ones,  and  the  covariances  become 
smoother  for  higher  altitudes  because  of  the  upward  continuation  effect. 
These  two  statements  might  sound  trivial  but  they  are  important, 
and  represent  guidelines  for  the  choice  of  grid  distances  As  and  At> 

The  reason  is  the  following:  Approximate  covariances  are  obtained  by 
a differentiation  - interpolation  procedure  applied  to  some  interpolating 
function  uniquely  defined  by  the  covariances  at  the  grid  points.  If  the 
interpolating  function  is  an  element  of  Kj  [d],  the  space  of  continuous 
functions  with  domain  D,  having  quadratically  integrable  first  derivatives, 
it  can  be  shown  that  the  maximum  approximation  error  depends  on  the 
i + l)th  power  of  the  grid  spacing  multiplied  by  the  maximum  absolute 
value  of  the  (X  + l)th  order  derivative  of  the  original  function.  Each 
differentiation,  if  it  exists,  decreases  the  dependence  on  the  grid  spacing 
by  one  power  (Sttnkel,  1978*.  p*  22),  (one-dimensional  case).  Consequently, 
the  grid  has  to  be  denser  for  small  spherical  distances  ^ than  for  larger 
ones  and  denser  for  small  altitudes  than  for  higher  ones. 
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Next  we  have  to  face  the  problem  of  which  kind  of  interpolation 
function  we  should  choose.  Five  requirements  have  to  be  fulfilled  by 
the  interpolation  function: 

a)  It  should  be  as  simple  as  possible  in  order  to  keep  the 
computation  time  at  a minimum. 

b)  It  should  be  as  smooth  as  possible  and  admit  as  many 
continous  derivatives  as  possible. 

c)  It  should  be  as  local  as  possible  in  order  to  guarantee  a 
good  fit  to  the  true  covariance  function  not  only  for  large, 
but  also  for  small  f -values. 

d)  It  should  disturb  the  true  spectrum  and  the  actual  degree 
variances,  as  little  as  possible. 

e)  It  should  satisfy  the  Laplac e-equation  with  reasonable  accuracy 
in  order  to  guarantee  harmonicity. 

Needless  So  say,  these  five  requirements  exclude  each  other: 
analytic  expressions  meet  the  requirements  (b)-(e),  but  fail  at  (a); 
single  higher  degree  polynomials  fulfill  partly  (a)  and  (b),  but  violate 
(c)-(e);  step  functions  and  piecewise  linear  functions  comply  with  (a) 
and  (c),  but  violate  (b),  (d)  and  (e).  Obviously  we  have  to  make  a 
compromise  between  competing  functions.  The  only  interpolation 
function  which  meets  all  five  requirements  sufficiently  well  is  the 
winner  of  the  game  - the  bicubic  spline  function:  it  is  relatively 
simple,  is  twice  continuously  differentiable  with  respect  to  each 
independent  variable,  has  very  local  features  although  it  doesn't  have  a 
compact  support;  it  reproduces  uie  actual  degree  variances  with  only 
little  noise  (Sunkel,  1978a,  p.  29a),  and  it  also  allows  the  Lapiace- 
equation  to  be  fulfilled  with  reasonable  accuracy  provided  the  grid  is 
sufficiently  dense. 


If  we  were  interested  in  covariances  involving  not  more  than 
three  differentiations  of  the  covariance  function  with  respect  to  s and  t, 
a single  bicubic  spline  representation  of  the  disturbing  potential's 
spatial  covariance  function  would  be  sufficient/  ' However,  we  intend  to 
derive  approximations  of  covariance  expressions  up  to  second  order 
derivatives  (at  both  points  P and  Q)  and  this  corresponds  to  fourth  order 
derivatives  with  respect  to  s or  t.  Consequently,  a single  spline 
representation  of  cov  (T,  T)  is  not  sufficient  anymore,  and  we  need  at 
least  two  more  representations.  We  have  chosen  cov  (Dr  T,  Dr,  T)  and 
cov  (Dt  T,  Dt  T)  = cov  (D*  T,  T)  in  order  to  guarantee  at  least  linear 
(and  at  most  cubic)  interpolation  within  the  grid.  These  three  spline 
representations  are  necessary  and  sufficient  for  our  purposes.  Since 
the  gravity  anomaly  is  the  most  frequently  used  quantity  in  physical 
geodesy,  we  have  extended  the  covariance  representations  and  included  also 
the  covariance  representation  of  the  gravity  anomaly  covariance  function. 

It  should,  however,  be  pointed  out  that  this  is  not  necessary  and  has  been 
done  only  for  reasons  of  convenience. 

The  four  covariance  functions,  cov  (T,  T),  cov  (Dr  T,  Dr,  T), 
cov  (D*  T,  T)  and  cov  (Ag.  Ag)  depend  only  on  s and  t and  so  do  their 
bicubic  spline  representations.  All  other  covariance  expressions  up  to 
(at  least)  second  order  derivatives  of  the  disturbing  potential  can  be 
obtained  from  linear  operations  applied  to  these  four  covariance  representa- 
tions. This  is  the  very  philosophy  behind  covariance  approximations  which 
helps  us  to  reduce  drastically  the  CPU -time. 


(*)...  This  follows  from  the  fact  that  the  3rd  derivative  of  a cubic  spline  is 
a step  function. 
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3.  Covariance  Grid  Point  Value* 

Since  a bicubic  spline  representation  of  the  covariance  function 
will  be  used,  it  is  necessary  to  determine  all  derivatives  of  the  covariance 
function  which,  together  with  their  function  values  at  all  grid  points,  make 
the  spline  representation  unique.  Let  D(s,  t)  be  a rectangular  domain 
divided  into  I*J  subrectangles  by  the  grid  points 

•o  * »t  * «i  . *<>  * t,  a t, 

with  (s,,  t|)  grid  point  coordinates.  Then  a bicubic  spline  on  D is 
uniquely  defined  by  (Sunkel,  1978a,  p.  74) 

*)  l * 0,  ...,I|  j = 0,  ...,J  ...function  values  at  the  grid  points, 

b)  D#ffl  * * * 0*  II  j * 0, ...,  J ...first  derivatives  of  the  function  f 

with  respect  to  s, 

c)  • * * 0,  ...,I|  j * 0,  J ...first  derivatives  of  the  function  f 

with  respect  to  t, 

d)  D|t{tj  , i = 0,  I)  j * 0,  J ...  second  derivatives  of  the  function  f 

with  respect  to  st. 

Here,  f stands  for  one  of  the  four  covariance  functions  cov(T,  T), 
cov(Df  T,  D,,  T),  cov  (D*  T,  T),  cov(Ag,  Ag).  What  we  need  are  expressions 
for  first  order  partial  derivatives  of  all  4 covariances  with  respect  to 
their  independent  variables  s,  t.  Since  the  subroutine  COVAX  of  C.C. 

T scheming  can  provide  us  with  values  for  these  quantities  implicitly,  we 
intend  to  express  the  derivatives  in  terms  of  COVAX  covariance  outputs. 

In  the  next  sections  we  derive  these  formulas. 

3.  1 Derivatives  of  cov(T,  T)  with  respect  to  s and  t 

We  recall  the  definition  (2 -1)  of  the  disturbing  potential  covariance 
function  which  we  write  in  a simplified  notation  as 
cov(T.T)  = Ek*  s«+*  P,(t), 


the  partial  derivative  with  respect  to  s is  trivial: 

D,  cov(T,  T)  = Z K (»  + !)•*  P»(*)- 

This  expression  can  be  easily  related  to  cov(-”r"Dr  T,  T) 
by  considering 

cov(DrT,T)  = -i-Sk,(n+  l)s“+iPB(t), 
and  we  obtain 

Df  cov(T,  T)  = cov  ( --L_DrT,  T).  (3-1) 

Even  simpler  is  the  calculation  of  Dt  cov(T,  T)J 

Dtcov(T,T)  = Sk,  s *+1  DtPn  (t).  (3-2) 

(All  closed  expressions  which  we  don't  write  down  explicitly  here,  can  be 
found  in  Tscherning  (1976).  Their  calculations  are  part  of  the  subroutine 
COVAX. ) 

Similarly  we  obtain 

D.  Dt  cov(T,  T)  = -jp  cov(-  -L  D„  T,  Dt  T).  (3-3) 

3.2  Derivatives  of  cov  (Dr  T,  D.,  T)  with  respect  to  s and  t 
By  straightforward  derivation  we  find 
cov(Dr  T,  D,.T)  = jL,  sk,  (n+1)2  s»+1Pa(t). 

Recalling  the  definition  of  s,  we  can  also  write 

cov(Dr  T,  Dr.T)  = -§—  Sk,  (n+1)2  s*+JP,(t). 

Applying  the  differentiation  operator  D(  we  obtain 

D.  cov(Dr  T,  Dr,T)  = -A*  Eh.  (n+1)1  (n+2)s  •+»  P,  (t), 
which  can  easily  be  shown  to  be  equivalent  to 

D,  cov(D,  T,  Dr,  T)  = (h^  cov  A-  D,  T.  D$,  t).  (3-4) 
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The  differentiation  operator  Dt  applied  to  cov(Dr  T,  Dr,  T)  yield* 

Dt  cov(DrT,  Dr,  T)  = • E (n+1)1  s *+l  D%  P,  (t)  (3-5) 

and  similarly  we  obtain  applying  Dt  to  equ.  (3-4) 

D#Dt  cov(DfT,  D,,  T)  = cov  1_  Dr  Dt  T,  D*,  T^.  (3-6) 

3.3  Derivatives  of  cov(D*  T,  T)  with  reject  to  s and  t 

All  covariance  derivatives  of  this  type  are  closely  related  to  equation 


(3  — 1)|  we  obtain 

cov(D*  T,  T)  = Ek,  s ■ +,D*P#  (t),  (3-7) 

D,  cov(D*  T,  T)  = cov(  -Id,T,  D^T),  (3-8) 

Dt  cov(D*  T,  T)  = Ek|  ( ,+1D*  P,(t),  (3-9) 

D,Dt  cov(D*  T,  T)  = £i  cov  ( - D,DtT,  D*  T).  (3-10) 


3. 4 Derivatives  of  cov(Ag.  Ag)  with  respect  to  s and  t 

All  covariance  derivatives  of  this  type  can  be  written  down  immediately 
replacing  T by  Ag  in  section  3.1: 

D,cov(Ag.  Ag)  = - JLcov(Dr  Ag.  Ag).  (3-11) 

Dt  cov(Ag,  Ag)  = cov(Ag.  Dt  Ag).  (3-12) 

D,D%  cov(Ag,  Ag)  1 - f cov(D,Ag.  DtAg).  (3-13) 

All  covariances  given  above  can  be  calculated  by  a simplified  version 
of  COVAX.  The  modifications  are  described  in  the  listing  of  the  program 
COVNET. 
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4.  Covariance  Expressions 


All  covariance  expressions  which  can  be  derived  from  the  bicubic 
spline  approximation  model  are  linear  combinations  of  the  grid  point  covar- 
iances listed  in  chapter  3.  For  example,  if  the  covariance  cov(T*,  Tt)is 
desired.it  is  necessary  to  calculate  the  corresponding  s and  t 
values  and  find  the  element  indices  (i,j)  such  that  s,  % s 2 s|+l  and 
t,  >t  itJ+1.  On  each  rectangle  [s(l  s1+1i  t](  tJ+1]  the  spatial  covariance 
function  is  approximated  by  a bicubic  polynomial  with  16  coefficients 
being  linearly  related  to  cov(T,  T),  D,  cov(T,  T),  D*  cov(T,  T)  and  D,Dtcov(T,  T) 
at  the  four  comers  of  the  corresponding  rectangle.  To  find  cov(T,  T)  for 
s and  t values  not  coinciding  with  one  of  the  comer  values,  is  a simple 
problem  of  interpolation.  (To  find  cov(T,  D%  T)  for  the  same  point  (s,  t) 
would  be  a differentiation  and  interpolation  problem. ) 

The  general  form  of  the  bicubic  polynomial  defined  on  the  rectangular 
domain  la,,  s,+l;  t,,  t,  + J is 

to)  j ■>  (ij)  , 

f (•.  t)  = £ t s i (•-•t)*  (4-1) 

k=0  i =0 

The  16  coefficients  result  from  linear  transformations  of  f,  D(f,  Dtf  and 
DtDtf  taken  at  the  4 comers  of  the  rectangle: 

ft)  * * o’  F H 


with 


G = 1 G*  G.  = "3/«2  Z,K  G = 3/*2  “2/*S 

0 G2  ’ Gl  -2/g  1/g*  • G*  -1/g  1/g*  » 


...  -3/gi  2/gJ 


10  A _ 0 0 

0 lj*  0 " [o  oj  * g!  “ ■*+*“■« 


(4-2) 


G 8 = tJ+i  _ J 


Fi+i,i  Ft+t,j+i 


*•1  Dt 

D.Dtfu 


Simple  matrix  multiplications  yield 


I F»H|  + F,Hj 

Gl  F,  +Gj  F,  [*  (Gf  F | + Gj  F,)  H, 


• +(GlFs  + Gj  F4)H,  . 


(4-3) 


Before  writing  down  the  necessary  derivatives  of  the  bicubic 
polynomial,  (4-1)  is  simplified  by  introducing  reduced  variables 


s and  t. 


s : = s - s. 


t s = t - t, 


so  that  (4-1)  can  be  abbreviated  (suppressing  the  element  indices  i,  j) 
and  has  the  form 


f(M)  v **• 


(4-1)' 


A general  derivative  of  (4-1)'  can  be  given  by  the  formula 


D*f(l,  I)  * 


* t* 


•»••»!  b £ /ky^\i 

ww 


(4-4) 


with  the  multi -index  «<  = («,,  ot).  l#l  = On  + ft.  0 * «„  <*  a 3, 

This  nice  expression,  however,  is  far  from  optimal  with  regard  to 
CPU-time  savings.  Therefore,  all  partial  derivatives  of  f actually 
used  in  the  computations  (indicated  by  subscripts)  follow: 


f(I,  i)  * aM  + i (a0|  + l (aM  + \ a0,)) 

+ * ((»»«  + f(»u  ♦ f(au  + Sa„))) 

+ I ((»ao  + K»ai  + *(*ii  + fa,,))) 

4 1 (aM  + f(a„  U(aB  + fa,,))))). 


(4-5a) 
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f,(l.  1)  3 *10  + * (*u  + i (*u  + **u)) 

+ i (2(a,0  + f (»jt  + t (a„  + ta„))) 

+ i‘3  (a,0  + ^ (*si  + ^ (a„  + ^*js))))» 

(i.  i)  = 2 (*jo  + i (»ai  + M*u  + *a„)) 

+ i*3  (a,0  + t (a„  + l (a„  + ta„)))), 

f%(I,  l)  = a0,  + i (a ||  + s (a, | + ia„)) 

+ t (2  (aM  + a (a„  + S (a„  + B a,,))) 

+ t*3  (a0,  + i (an  + I (a15  + 1 a, ,)))), 

U%  (*.  t)  * 2 (aM  + i (au  + B (a„  + la,,)) 
+ 1*3  (a0,  + 1 (a |,  + B (a„  + I a„)))), 

f.%  (5.  t)  * »»»  + t (2  au  + t *3  iu) 

+ I (2  (a, | + t (2  a„  + t-3  a,,)) 

+ i-3  (a, | + l (2  a„  + i-3  a,,))), 


(4-5b) 


(4 -5c) 


(4-5d) 


(4-5e) 


(4-5f) 


f(tt  (•»  t)  = 2 (a,|  + T (2  a„  + f*3  a„) 
+ 1*3  (a„  + i (2  a„  + t*3  a,,))). 


(4-5g) 


The  reader  ia  invited  to  verify  that  the  function  values  together 
with  the  partial  derivatives  fa , f%  and  fa,  are  reproduced  at  the  four 
corners  of  the  rectangle  (hint!  use  (4-5),  (4-2)  and  (4-3). 

In  all  covariance  expressions  derivatives  are  taken  with  respect  to 
the  spherical  coordinates  (r,  p,  X),  (r\  p \ \ *)  of  the  points  P and  Q 
and  not  only  with  respect  to  s and  t.  According  to  the  chain  rule  for 
differentiation  we  need  to  know  the  partial  derivatives  of  s with  respect 
to  r,  r'  and  of  t with  respect  to  p,  p \ x'i  the  latter  are  listed 


explicitly  in  (T scheming,  1976,  pp.  18,  19).  To  find  the  redial 
derivatives  of  s using  (2 -2b)  is  also  trivial.  The  result  is 


Furthermore,  whenever  calculating  derivatives  of  f with  respect  to 
latitude  and/or  longitude,  it  is  necessary  to  be  aware  of  thefact  that  an  f€tf  -th 
order  derivative  consists  of  a linear  combination  of  all  derivatives  of  f 
with  respect  to  t up  to  the  lot  -th  order. 


with  mt,  i * 1 n any  one  of  fo,  X,«p\ x'|  and  * = <*„ ). 

The  coefficients  C|,  i 1 1 consist  of  derivatives  of  t with  respect 

to  fp,  1,  p x'|.  General  formulas  are  easy  to  derive  and  can  be  found 
in  (T  sc  homing,  1976,  p.  16). 

On  the  neat  pages  we  give  all  basic  covariance  expressions  in  their 
approximated  version  are  given.  The  following  abbreviations  will  be  used: 

fl.  . . bicubic  spline  representation  of  cov  (T,  T), 

f2. . . bicubic  spline  representation  of  cov  (Dr  T,  DK  T), 

£3...  bicubic  spline  representation  of  cov  (D*  T,  T), 

f4. . . bicubic  spline  representation  of  cov  (Ag.  Ag). 


Several  times  the  subsequent  identities  were  used; 


cov 


cov 


After  having  established  the  relations  between  basic  quantities,  the 
"long  march"  of  deriving  all  covariance  expressions  begins.  The  results 
are  listed  below. 


cov  (T,  T)  = fl 


cov  (T,  D‘ 


cov  (T,  D 


cov 


I 
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cov  (DrT,  Dr,  T)  = £2 
cov  (Dr  T,  Ag)  = -£2  + 2 £1,  tt' 

cov  (Dr  T,  D,.  Ag)  ■ (£2.«-  2 £a  - 2 £1,  —,)  “ 

cov  (DrT.  D^,  T)  = -£2#  -Jr 

(4-10) 

cov  (DrT.  D T)  = -£lit  *J"  c, 

l 

cov  (Df  T,  Dr,  T)  = £2t  c, 

cov  (Dr T,  Ag)  = (-£2t  + 2 £l-%  -Jp#)  c» 

cov  (Df  T,  T)  = - (£3,  c,  + £1>%  C|)  4 

cov  (Ag.  Ag)  = £4 
cov  (Ag.  Dr.  Ag)  - - £4,  -Jr 
cov  (Ag,  D$.  T)  = cov  (DrT,  D,.  Ag) 

cov  (Ag.  T)  * (£li%  • - 2 £1%)  -p-  c,  (4-11) 

cov  (Ag.  Dr,  T)  = (-£2t  + 2 £1#%  -pp  c, 

cov  (Ag.  Ag)  = £4t  c, 

cov  (Ag.  T)  = [ (£3.  • - 2 f3)  c,  + (fl,ta  - 2 fl%)  ct]  -j- 

cov  (Df  Ag.  Dr,  Ag)  = (£4#<a  + £4, ) 4 

cov  (Dr  Ag.  Dj,  T)  = (£2i#.  - £2,  - (f2  -J-+  fl.  “))£■, 

cov  (Df  Ag.  D T)  = (-r'  f2%  + fl,t-J-  - -§-£lt  ) "T  c» 

, . 1 (*-12> 
cov  (Dr  Ag.  D^Dr,T)  * (£2#%a  - 2 £1  #l?p-  2 £2t)  — c, 

cov  (D,  Ag.  Ag)  = -£4it-J-  ca 

cov  (D,  Ag.  D#i  D#jT)  = [ (2  £3  - £3..  a1)  c,  + (-rr'  £2t  + £l,,a  + 2£l^)ct]  jr 


[ 


cov  (D| 
cov  (D‘ 
cov  (D‘ 
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cov  (D; 
cov  (D: 


T.  Dj.  T)  = (f2gga  + f2g)  . 

T.  D.  T)  = (r'  £2t  +~  fl% ) c» 

i 

T,  ^ Df/  T)  = - f2g  t -p~  C|  (4-13) 

T.  D#i  Ag)  = (f2gts  - 2 f2t  - 2 flg%  £,) -Fcl 

T*  Dtf  D T)  = [ (f3..  • + 2 £3,)  • c,  + (rr'  £2t  + £1,  %a)  C|)  ~t 

* "i 


cov  (Dw  T,  Dv  T)  = £3  c,  + flt  c, 

cov  (D„1  T.  DaiD,,T)  = -<f3.  c,  + f»„  c.)  -7  (4.M) 

cov  (D^1  T,  D0j  As)  - [ (f 3,  s - 2 f3)  c,  + (£1, ,«  - 2 £lt)  ct]  TV 

cov  T,  D#i  Daj  T)  = £3t  c,  + £3  c,  + flt  cx 

cov  (Dai  Dr  T,  Dr,  T)  = f2tt  c,  + £2t  cx 

cov  (Djj^  DrT,  D^  Ag)  = -(£2tt  - 2 £3g  Tr*)  c*  " ^ ^1«»  rTF*)  ci  (4-15) 

cov  (D#i  Dr  T,  D#j  T)  = (£3gt  c,  + f3g  c,  + flgt  c,)  ar 


cov  (D#  Ag.  D#  Ag)  = *4U  c,  + f4t  cx 

(4-16) 

cov  (DWi  Ag.  Dtfj  D#j  T)  = [ f3lt  cs  + f3g  c,  +£ig%c,)a  -2  (£3%c,  + £3  c,  + fltc,)]  — 

cov  (D#i  DtfjT.  Dttj  D*4  T)  = f3lt  €4  + £3t  c,  + £3  c2  + flt  c,  (4-17) 

All  coefficients  c,,  i = 1, . . . , 4 occurring  in  the  foregoing 
exprea8iona  are  identical  with  the  coefficienta  of  Kt  (the  partial  derivativea 
of  the  covariance  function  with  reapect  to  t)  explicitly  liated  in  (Tacheming, 

1976,  pp.  17,  18,  19).  The  equationa  derived  in  thia  chapter  permit  the 
computation  of  all  covariance  expreaaiona  for  aecond  and  lower  order 
derivativea  of  the  diaturbing  potential. 
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5.  Program  Guid*li«« 

The  very  idea  behind  covariance  approximations,  as  stated  in  the 
introduction,  can  be  sketched  as  follows: 

a)  generates  sufficiently  dense  network  of  basic  covariances  (here: 
cov  (T,T).  cov  (DtT,  Dr,T),  cov  (D*  T,  T),  cov  (Ag.  Ag)  at  the  grid 
point,  of  a regular  rectangular  grid);  calculate  the  coefficients  of 
the  interpolating  functionand  store  this  information  once  on  some 

file. 

b)  compute  all  other  covariances  by  a simple  differentiation  - interpolation 
procedure  applied  to  the  interpolation  function  chosen  in  (a). 

The  program  is  organised  accordingly:  it  consists  basically  of  two  parts, 
the  covariance  network  generation  program  COVNET  and  the  covariance 
approximation  program  COVAPP. 

5*  * The  covariance  network  generation  program  COVNET 

According  to  its  purpose  the  program  performs  successively  the 
following  tasks: 

a)  it  reads  the  parameters  of  the  covariance  function  which  coincide 
with  the  parameters  used  in  Tscherning's  COVAX  subroutine. 

b)  It  r.dl  the  pftrftmeters  of  th.  cov.ri.nc.  grid  .nd  c.lcul»t«.  the 
gnd  v»lu«,.  Th.  folioring  figure  m»y  help  to  m»k.  th.  meening 
of  the  grid  parameters  more  clear: 


I 
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The  grid  in  s direction  ie  equi spaced,  though  not  with  respect  to 
the  corresponding  height.  This  is  a consequence  of  the  mapping 
• = R^/(rr#).  However,  since  the  values  of  s differ  only  slightly 
from  1 for  points  in  the  close  neighborhood  of  the  'irth,  the  s - grid 
with  constant  grid  spacing  A>  corresponds  to  a grid-height  with 
slightly  increasing  grid  spacing  Ah  when  h increases.  It  would  have 
been  also  possible  to  choose  an  s - grid  such  that  the  grid-height 
spacing  Ah  is  kept  constant;  however,  this  would  have  complicated 
the  retrieval  of  the  element  indices  corresponding  to  specified  (s,  t)  - 
values  and,  moreover,  is  of  no  advantage  at  all. 

After  the  grid  calculation,  the  program  performs  the 


c)  calculation  of  the  covariances  at  the  grid  points.  This  is  done  with 
a simplified  version  of  COVAX,  called  COVAXS.  The  necessary 
equations  are  given  in  chapter  3.  It  was  decided  to  calculate  the 
covariances  at  all  grid  points  and  not  only  those  which  define  the 
bicubic  spline  uniquely.  Thereby  obtaining  numerical  stabilization 
and  some  accuracy  improvement.  Finally  the  polynomial  coefficients 
are  calculated  for  all  elements  using  the  subroutine  BILDEC. 


5.2  The  covariance  approximation  program  COVAPP 


This  program  renders  in  principle  the  same  as  COVAX  does;  it 
provides  covariances  between  14  different  quantities  which  can  be 
described  as  linear  functionals  applied  on  T.  These  quantities  are 
described  in  (Tscheming,  1976);  however,  for  the  sake  of  completeness 


1 


height  anomaly 


2 

3 

4 

5 

6 

7 

8 

9 

10 
LI 
L2 
L3 
14 


— D T 

r r 


Ag  = -Dr  T-  -T 
-Dr  Ag 


d;  t 

l 


5 = - rZ  Dcp  T 


-7 


'yCO  Sep 


r"  D<pDr  T 


D\T 


1 


DXDrT 


r cos  cp 

“Y~  DcP  Ag 


:DX  Ag 


D 2 T 


r cos  op 
1 

r2-  "cp 
^'cos^'  D<PD*  T 


1 


CO 


-Di T 


9p 


radial  derivative  of  T divided  by  r 
free  air  gravity  anomaly 

radial  component  of  the  gravity  anomaly  gradient 
rr  - component  of  the  second  order  tensor  of  T 
N-S  component  of  the  vertical  deflection 
E-W  component  of  the  vertical  deflection 
rep  - component  of  the  second  order  tensor  of  T 

r\  - component  of  the  second  order  tensor  of  T 

N-S  component  of  the  gravity  anomaly  gradient 
E-W  component  of  the  gravity  anomaly  gradient 
epep  - component  of  the  second  order  tensor  of  T 

cpX  - component  of  the  second  order  tensor  of  T 

XX  - component  of  the  second  order  tensor  of  T 


Note  that  the  codes  of  (8,  9)  and  (10,  11)  have  been  interchanged  compared 
with  (T scheming,  1976,  pp.  2,  3).  This  change  simplified  the  program's 
logic.  The  codes  of  the  linear  functionals  on  T taken  at  the  points  P and  Q 
are  transferred  to  COVAPP  by  KR1  and  KR2.  The  coordinate  informations 
about  P and  Q are  transferred  by  the  vector  CR.  The  (s,  t)  coordinates 
are  checked  for  their  position  within  the  admissible  ranges  specified 
by  the  grid,  an  error  message  is  returned  if  s and/or  t is  out  of  the 
range.  Before  the  actual  calculation  of  the  covariance  COV,  all 
necessary  polynomial  coefficients  corresponding  to  the  element  in 
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consideration  are  made  available  by  the  subroutine  GET.  The  partial 
derivatives  of  the  bicubic  polynomial* s)  combined  with  interpolation  is 
rendered  by  the  function  BSFC.  The  calculated  covariance  COV  is 
returned  to  the  calling  program. 

In  order  to  provide  the  user  with  a fast  comparison  between  COVAPP 
and  COVAX,  we  have  designed  a program  called  TEST.  Detailed 
descriptions  are  integrate  parts  of  each  program. 


6.  Conclusions 

In  this  report  an  approximation  model  of  the  spatial  covariance  function 
of  the  disturbing  potential  ha.  been  described.  The  idea  was  to  generate 
once  a basic  covariance  network  which  is  chosen  sufficiently  dense  such  that  it 
permits  the  calculation  of  all  covariances  (up  to  and  including  second- 
order  derivatives)  by  applying  a simple  differentiation  - interpolation 
procedure  to  an  appropriate  interpolation  function.  It  turned  out 
that  the  best  possible  function  (if  "best"  is  defined  by  criteria  like 
simplicity,  smoothness,  localness,  etc.)  for  such  a purpose  is  a bicubic 
spline  function.  The  procedure  described  here  is  able  to  provide  the  same 
covariances  as  the  subroutine  COVAX  (Tscheming,  1976),  with 
arbitrary  high  accuracy  (depending  on  the  grid  spacing). 

The  goal  was  to  cut  down  significantly  the  computation  time  of  covariance 
expressions.  Covariance  expressions  for  second  and  lower  order 
derivative,  of  the  disturbing  potential  are  available  from  COVAPP,  a 
covariance  approximation  program  at  a CPU-time  level  between  1.8  10*«  sec 
(cov  (T,  T))  and  5.2  10*4  sec  (second  order  derivatives)  calculated  on  an 
IBM  370  system.  This  corresponds  to  a tenfold  gain  in  calculation  speed 
compared  with  COVAX.  which  w.  believe  to  be  a significant  improvement. 
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C c 
C C 0 V If  E T COMPUTES  THE  BASIC  NETWORK  OF  COVARIANCES  USING  A C 
C SIMPLIFIED  VEHSION  OF  C.C.  TSCHEHN I NGS  COVAX.  CAIXEDC 
C COVAXS,  AND  CALCULATES  THE  POLYNOMIAL  COEFFICIENTS  C 
C FOR  ALL  ELEMENTS  OF  THE  GRID.  C 
C C 
C THE  PROCRAM  IS  DESIGNED  FOR  VARIOUS  SOURCES  OF  GEODETIC  DATA  i C 
C SURFACE,  AIRBORNE  AND  SATELLITE  DATA.  THE  DIFFERENT  ALTITUDES  C 
C ABOVE  THE  EARTH  TO  WHICH  THE  DATA  REFER  MAKE  IT  NECESSARY  TO  C 
C ADOPT  A NON-UNIFORM  GRID.  THE  GRID  USED  FOR  HANDLING  SURFACE  C 
C AND  AIRBORNE  DATA  ONLY  WILL  BE  VERY  DENSE  ( ESPEC I ALLY  FOR  SMALL  C 
C SPHERICAL  DISTANCES)  COMPARED  WITH  THE  CRID  USED  FOR  SATELLITE/  C 


C SURFACE  (AIRBORNE)  DATA  OR  FOR  SATELLITE/SATELLITE  DATA.  IN  ORDERC 
C TO  BE  APPLICABLE  FOR  ALL  DIFFERENT  KINDS  OF  DATA  SIMULTANEOUSLY,  C 
C THE  PROGRAM  HAS  BEEN  DESIGNED  IN  SUCH  A WAY  THAT  THE  USER  CAN  C 
C DEFINE  0 RANGES  IN  HEIGHT  AND  5 RANGES  IN  SPHERICAL.  DISTANCE.  C 
C EACH  WITH  VARIABLE  GRID  SPACING.  C 


C C 

rrnrf.rrrrrrm,(Tf:r.rrrrr.rr(’rrrnrrrriwrrrrrrrrrr<'iwf'f'i'r'rrr('rprrri'<'ri'i'<vi' 


**********  INPUT  ********** 

COVARIANCE  MODEL  PARAMETERS: 


C 

C 

C 

C 

C 


C 

RBRE3  . . . SQUARE  OF  THE  RATIO  BETWEEN  THE  RADIUS  OF  C 

THE  BJERHANMAR  SPHERE  AND  THE  MEAN  RADIUS  C 
OF  THE  EARTH  (RE)  C 

RE  ...  MEAN  RADIUS  OF  THE  EARTH  C 

A ...  SCALING  FACTOR  OF  THE  GRAVITY  ANOMALY  DE-  C 

GREEK  VARIANCE  NODEL,  UNIT  ...  ( K/8EC)**4  C 
KK3).  KI< 4)  ...  INTEGERS  K2  AND  KB  IN  FORMULA  (17)  IN  C 

TSCHERNING  ( 1976)  C 

KI(8)  ...  NUMBER  OF  DEGREE  VARIANCE  NODEL  AS  SPEC I-  C 

FIED  IN  TSCHF.RN I NG  ( 1976)  i POSSIBLE  NUN-  C 
BERS  t 1,  2,  3 C 

N ...  DECREE  VARIANCES  UP  TO  AND  INCLUDING  DE-  C 

GM’F.E  N ARE  EITHER  SET  EQUAL  TO  ZERO  OR  C 
REPLACED  BY  EMPIRICAL  DECREE  VARIANCES.  C 

IN  THE  FIRST  CASE  THE  LOGICAL  VARIABLE  C 

LOCAL  MUST  BE  .TRUE.,  IN  THE  SECOND  CASE  C 
.FALSE.  C 

LOCAL  ...  LOGICAL  VARIABLE  AS  SPECIFIED  ABOVE  C 

SI OMAN ( I) , I* l.N+t  ...  EMPIRICAL  ANOMALY  DECREE  VARIANCES,  UNITS  C 

OF  NGAL**2  REPLACING  THE  NODEL  DECREE  C 

VARIANCES  UP  TO  AND  INCLUDING  DECREE  N.  C 
( S IGMAC  K+ 1 ) CORRESPONDS  TO  THE  DECREE  C 

VARIANCE  OF  DECREE  K>  C 

C 

COVARIANCE  GRID  PARAMETERS  : C 


C(MH,MD,K. 16)  ... 


NDR  ... 

W>(  I)  , 1*1, NDR  ... 
Hl( I) , Is l, NDR  ... 
DH( I) , I* I , NDR  ... 

S«(  I) , I- I, NDR  . . . 

8 1 ( I)  , 1*1, NDR  . . . 
DS(  I)  , 1*1, NDR  . . . 

NDPS( I). I* I, NDR  ... 

DPSKI.J)  ... 


4-DIM.  ARRAY  OF  POLYNOMIAL  COEFFICIENTS  > C 
K«l  . . .COV(T.T)-COEFFIClENTS,  C 
K«2  . . .COV(  DR(  T) , DR’  ( T) ) -COEFFICIENTS,  C 
K-3  . . COV(DT(T>  ,DT(T»> -COEFFICIENTS,  C 
K-4  ...COV( DELTA  G, DELTA  C) -COEFFICIENTS  C 

C 

NUMBER  OF  RANGES  IN  RADIAL  DIRECTION  C 
( NDR. GE.  I. AND. NDR. IJS. 3)  C 
LOWER  LIMIT  OF  HEIGHT  RANGE  I C 
UPPER  LIMtT  OF  HEIGHT  RANGE  I C 
STARTING  CRID  SPACING  IN  HEIGHT,  RANCE  I C 

C 

LOWER  LIMITS  OF  RADIAL  RANGES  IN  8 C 
UPPER  LIMITS  OF  RADIAL  RANGES  IN  S C 
GRID  DISTANCES  IN  S-DIRECTION  C 

C 

NUMBER  OF  RANGES  IN  8PHER I CAL  DISTANCE  C 
( NDP8( I) .GE. I . AND. NDPS(  I ) . LE . 8 > C 
SPHERICAL  DISTANCE  CRID  SPACINCS  (I  STAN DSC 
FOR  HEIGHT  RANGES,  J FOR  SPH.  DIST.  RANCESC 


M 

! 

i 


( 


I 
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ndpshi.j)  ... 

P8IS(I,J)  ... 
P8IKI)  ... 

P8< I , K)  ... 

IACC(I),  JACC(I.J) 


NUMBER  or  SPHERICAL  DISTANCE  CRID  SPACING8C 
( 1 AND  J AS  ABOVE)  C 

SPHERICAL  DISTANCE  RANGE  STARTING  VALUES  C 
SPHERICAL  DI8TANCE  UPPER  LIHIT  VALUES  C 
EACH  ROW  VECTOR  CONTAINS  THE  COSINE  Of  THEC 
SPHERICAL  DISTANCE  OP  THE  CRID  POINTS  C 
INTEGER  VECTORS  USED  TO  AS8ICN  NUMBERS  TO  C 
GRID  POINTS  C 

C 

OUTPUT  »***»**#*»  C 


RV(  L)  . , 
DV<  I.J) 


VECTOR  Or  GRID  VALUES  IN  8 DIRECTION  C 
EACH  ROW  VECTOR  CONTAINS  THE  COSINE  OP  C 
THE  SPHERICAL  DISTANCE  OF  THE  GRID  POINTS  C 

C 


THE  PROCRAH  USES  THE  SUBROUTINES  CO VANS,  BILDEC  AND  BLOCXDATA.  C 

C 

C 

tr  COVNET  AND  TEST  ABE  RUN  SEPARATELY,  THE  STATEMENT  C 

• CALL  TEST*  C 

HAS  TO  BE  REPLACED  BY  A VR I TE- STATEMENT  SUCH.  THAT  ALL  QUANT I-  C 
TIES  ENCOUNTER  IN  THE  LABELED  COMMON  BLOCKS  /CHID/  AND  ✓CRIPAR'  C 
ARE  STORED  ON  PILE  I.  C 

C 

IN  CASE  YOU  LIKE  TO  CHOOSE  M INSTEAD  OP  8 AS  THE  NUMBER  OF  RADIALC 
RANGES  AND  N INSTEAD  OP  8 AS  THE  NUMUCR  OF  SPHERICAL  DISTANCE  C 
RANGES,  YOU  HAVE  TO  REPIJtCE  ALL  3 BY  M,  ALL  8 BY  N AND  ALL  6 BY  C 
NM  WITHIN  THE  LABELED  COMMON  BLOCK  /CRIPAR/.  THIS  CHANCE  HAS  TO  C 
BE  MADE  IN  COVNET,  TEST  AND  COVAPP.  C 

C 

IN  CASE  YOU  LIKE  TO  CHANGE  THE  DIMENSION  OF  THE  CRID  ARRAY  PROM  C 
COO, 83, 4,  16)  TO  C(KI,K2,4, 16)  , YOU  HAVE  TO  MAKE  THE  SAW  CORREC-C 
TIONS  ALSO  IN  TEST.  BILDEC  AND  GET.  BESIDES  THAT  THE  DDENSIONS  C 
OF  RVOO)  AND  DVO.83)  HAVE  TO  BE  NOW  RV(KI),  DVO.K2).  CHANGE  C 
THE  DIMENSION  OP  DV  ALSO  IN  TEST  AND  OP  THE  CORRESPONDING  ARRAY  C 
PS  IN  COVAPP.  CHARGE  ALSO  THE  VALUES  POR  HI  AND  HD  IN  THE  DATA  C 
STATEMENT  IN  COVNET.  C 


REFERENCES  « 


SUENKEL  < 1978) 


C 
© 

C 

C 

C 

C 1. 
C 

C a. 
C 

c 

c 

C 8. 

C 

C 

C 

C 

C 

c 

C 

c 

C 

c 

C 1. 


c 

c 8. 


C0VAX8  AND  COVAXN  DIPPER  PROM  C O V A X < TSCHER- 
NING,  1974)  BY  THE  FOLLOWING  MODIFICATIONS  i 

THE  COMPUTED  COVARIANCE  COV  IS  AN  ELEMENT  OP  THE  PARAMETER  LIST. 

THE  ENTRIES  COV AX.  COVBX  AND  COVCX  HAVE  BEEN  REPLACED  BY  VARIAB- 
LES SUCH.  THAT  (-I.COV)  CORRESPONDS  TO  THE  ENTRY  COVAX,  <3. COV) 
TO  THE  ENTRY  COVBX  AND  ( I. COV)  TO  COVCX. 

THE  CRITICAL  ALTITUDE  HHAX  IS  KEPT  CONSTANT  AND  EQUAL  TO  28  KM. 
IPIKKBXS  OR  I Kl( 8) BS  AND  Kl( 4) >KI(8) ) > AND  IF(KKO)  OR  Kl(4)> 
333) , LSUH  IS  .TRUE.,  OTHERWISE  .FALSE.  . 

IP  THE  FIRST  CONDITION  IS  NOT  FULFILLED,  L8UH  IS  . TRUE.  POR 
KHS)  OR  KH4)  > 233.  OTHERWISE  IT  IS  .FALSE.  . 

WHEN  LSUH  IS  .TRUE.,  WILL  THE  COVARIANCES  IN  ALTITUDES  HIGHER 
THAN  28  EH  BE  COMPUTED  BY  EVALUATING  THE  SUM  OP  THE  LEGENDRE  - 
SERIES  HAVING  N2-299  TERM, 

ADDITIONAL  8HVLIPICATI0N8  POR  C O V A X 8 t 

THE  INTEGER  ARRAYS  K19,  K21  AND  K28  HAVE  IDENTICALLY  THE  VALUE 
ZERO. 

THE  ARRAY  CM  HAS  IDENTICALLY  THE  VALUE  ONE. 

SINCE  COVAJB  MUST  ONLY  PROVIDE  DERIVATIVES  OP  THE  DISTURBING 
POTENTIAL  COVARIANCE  FUNCTION  IN  RADIAL  DIRECTION  AND  WITH  RE- 
SPECT TO  THE  COSINE  OP  THE  SPHERICAL  DISTANCE  BETWEEN  TWO  POINTS, 
THE  COMPUTATION  OP  THE  QUANTITIES  D<  I) D<34>  ( TOCHER- 


• • • 9 
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C NINC,  1976.  P.  57)  IS  NOT  NECESSARY.  FOR  THE  SANE  REASON  COVANS  C 

C NEEDS  ROT  TO  CONFUTE  THE  INTEGERS  SPECIFYING  THE  KINDS  OF  DIFFKR-C 

C ENT I AT I ON  AND  THE  DIFFERENTIATION  ITSELF  WITH  RESPECT  TO  THE  C 

C LATITUDES  AND  LONGITUDES.  THEREFORE.  THE  PART  BETWEEN  THE  STATE-  C 

C KENTS  #78  AND  #85  18  ROT  NECESSARY  EITHER.  STATEMENT  #85  HAS  TO  C 

C BE  REPLACED  BY  n 

C ‘85  COV*C!  NDI+ 1 ) * C 

C (TSCHERNING.  1976.  P.  68).  C 

C THIS  LIST  GIVES  ONLY  THE  HOST  E88ENTIAL  CHANGES  AND  PRO-  C 

C VIDES  THE  USER  A TOOL  TO  OUICKLY  ADJUST  COVAX  IN  ORDER  TO  BE  C 

C FULLY  APPLICABLE  FOR  COVNET  AND  TEST.  RESPECTIVELY.  C 

c C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCOCC 
IMPLICIT  REAL*8< A-H.O-Z) 

LOCI CAL  LOCAL 

COMMON  /CRID/C!36.56.4. 16)  /CRIPAR/RB2.  86(3),  DS(8>,  61(3),  DP8K 

* 3,5)  . PS180.6)  . P81 1(3)  , DV!S,56)  , RE.  IACC(3>.  JACCI! 3.5) 

* 

* /COM4/RV1 36) 

8ICMA6I366),  8ICNA(  366)  . KK85).  Nl. 


3,5).  PSISO.6),  PSI  M3)  , DV<3,5#>,  RE.  IACC(3),  JACCK ! 
. NDR  /COIQ/DUHNY!  2)  , DX.  DY.  I,  J 

/COH3/H6!3),  DH(  3) , Hl(  3) , JACC(3).  NDPS(3>  , NDP8K3.5) 


CR(  51) 


3661 

5662 

5663 

5664 

II 


6662 

6663 

6605 

6666 

6667 

6668 
6669 


6616 

6611 

6612 

6613 

6614 

6615 

6616 

6617 

6618 
6619 

C 

C 

C 


C 

C 

C 


* /CHCOV/CI!  12) 

* LOCAL 

DATA  RHOMBS. 437746 77 I D3/.NH.  ND/36.56/ 

FORMAT* 3D 1 5 . 1) 

FORMAT!  2D 15. 7. 4 15.  U» 

FORMAT!  12F6.2) 

FORMAT! 1615) 

FORMAT!  5(  DI6. 2, IS)) 

FORMAT!  I H .4X. 'COVARIANCE  FUNCTION  PARAMETERS  « * ,//,5X,  *1 
1DI5.7./.5X.  * A ' , IDI5.7./.5X,  'EM 3) ' , I3./.5X, 


- I3./.5X,  'KI(S)  •••••  | I0|F |OA|  Hi i • • • •••••! f 9 

FORMAT!  IH  .4X, 'EMPIRICAL  ANOMALY  DEGREE  VARIANCES  IN  UNITS  OF  HGAL 
***2  : ' .//.5X.25!  I2F6.2/) ) 

FORMAT! I H6.4X. 'CPU-TIME  NEEDED  FOR  THE  CALCULATION  OF  THE  COVARIAN 
»CE  GRID. . . ' .FI6.5/I 
FORMAT! IH  ,4X. 'LEVEL. . . ' , 13/) 

FORMAT!  IH  ,4X. 16D12.4) 

FORMAT!  IH  ,4X. 'NUMBER  OF  GRID  POINTS  IN  RADIAL  DIRECTION  EXCEEDS  T 
*HE  PRESCRIBED  LIMIT  IN  THE  DIMENSION  STATEMENT'/) 

FORMAT!  I H ,4X. ’NUMBER  OF  GRID  POINTS  IN  SPHERICAL  DISTANCE  EXCEEDS 
*THF.  PRESCRIBED  LIMIT  IN  THE  DIMENSION  STATEMENT'/) 

FORMAT!  I H6.4X, 'NUMBER  OF  RADIAL  RANGES  ! NDR)  . . . ’ . I3.//5X. * I 

* H6I 1)  HI! I)  DH!  I)  ' ,/.5X,  ' 

*  •/) 

FORMAT! IH  . 4X. 12, 2X, 3FI5 . I/) 

FORMAT!  I H6.4X. 'NUMBER  OF  SPHERICAL  DISTANCE  RANGES * //5X, • I ND 

*PS<  I)  ' ,/,5X,  ’ • ./) 

FORMAT!  IH  ,4X. I2.6X. 12) 

FORMAT!  1H6.4X, 'SPHERICAL  DISTANCE  GRID  VALUES  AND  NUMBER  OF  CORRES 
* POND INC  INTERVALS’ //8X, • I J DPSl(I.J)  NDP8K I , J) ’/) 

FORMAT!  IH  , 4X. 12, 13 . F 14. 6, 6X. 13) 

FORMAT!  1116, 4X, ’GRID  POINT  VALUES  IN  RADIAL  DIRECTION*/) 

FORMAT!  IH  .4X.BFI2.1) 

FORMAT!  I H6, 4X, ’GRID  POINT  VALUES  IN  SPHERICAL  DISTANCE’/) 

FORMAT!  I U6. 4 X. ’RADIAL  RANGE  NUMBER  . . . * . 12/) 

FORMAT!  IH  .4X.8FI6.6) 

RE*637 I . 6D3 

READ  AND  ECHO  COVARIANCE  PARAMETERS 

READ! 5, 5661)  RBRE2, A, Kl!3) ,KI(4) .XI (5) ,N, LOCAL 
WRITE!  6. 6661)  RBRE2,  A.  KI  ( 3)  . EH  4) . KU  5) . N 
IF! N.LT.2)  N-2 
Nl-N-M 

IF!  .NOT. LOCAL)  READ! 5, 8662)  (8IGMA6! I) , I-l.HI) 

IF! .NOT. LOCAL)  WRITE! 6,6662)  (8IGMA6!  I), !■ I.NI) 

READ  AND  ECHO  COVARIANCE  GRID  PARAMETERS 

READ! 8,5663)  NDR 

READ! 5.5666)  ! ! H6(  I ) , HI!  I ) . DH!  I) ) , I* 1 , NDR) 

READ! 5.3663)  ! NDPS!  I ) , I ■ I , NDR) 

DO  21  I ■ I , NDR 


, I3./.3X,  *N. 


13//) 


KK4) 


20 


SI 


C 

C 

C 


I 


a 


is 

19 


as 

3 

as 

as 


aa 


as 

c 

c 

c 


SDR* SOPS! 1) 

REAIKS.SSSS)  ((BMK  I.J).NDP8I!  I . J) ) . J*  I . SDR) 
WRITE!*. SSS9)  SDR 

WRITE!  A, SO IS)  ( ( I, S3!  I)  .HI!  I)  ,DS<  I)  > . !•  I .HSR) 
VRITt(  § i W 11) 

Sft]TE<«±6tia>  !!I,*DP8!I)>,I*1,*DR) 


CI!8)-A*RBa*l.SD-lS 
CI<  IS) -I 


CALCULATE  THE  ARRAYS  SS.DS.S1.RV. IA0C.P8I8.P8I l.DV.JAOC. JACCI 


I-S 

I AOS 
I*  Is  I 
I AC*  I AO  I 

ss(  n*iua/(as(  n+RE)/(BS(  d+re) 

RV(  IAC) *HS( D+RE 
I AC  1*1  ACM 

RV( IACI > *HS( I)sDH<  I) ♦RE 
8*RBa/RV(  I AC t ) /RV! IAC1) 

DS(  I) *8-88!  I) 

IAC- I AC  1 
8-8+08! I) 

R-RBa/8 

R*D8QRT!R) 

1F!R.GT.!H1<  I)+DH! I)+RE)>  GOTO  19 
IAC* IAC+ I 

IF( IAC.LE.NH)  GOTO  IB 
WRITE! 6 , 6SS7) 

STOP 

RV< IAC) *R 

goto  a 

Sl(  I)*8-D6!  I) 

IACC< I) ■ IAC 
IF! I. LT. SDR)  GOTO  I 
do  as  I* 1 , SDR 
P8I8< I, 1)*S.SDS 
DV<  I.D-l.SDS 
JAC*  I 

m>R*RDP8(  I) 

P8IL-S.SDS 

RlfR-S 

DO  as  J* I , SDH 
RDN-SDPSK  I, J) 

DPSK  I,  J)-DP8I(  I . J)/RHOH 

SRH-RHS+SDH 

JACCI(  I,J)*HHR 

DO  3 K* 1 , SDR 

JAC* JAC + I 

IF( JAC.LE.HD)  GOTO  2S 
WRITE!  S.6SS8) 

STOP 

P8IL*P8IL+DPSI!  I,J) 

DV(  I,  JAC)*DC08<P8IL> 

P8I8(  I.J*I)-P8IL 
JACC< I)* JAC 
P8IK  I)*PSIL 
WRITE!  *,*813) 

DO  aa  I* 1 ,RDR 
SDN*SDP8< I) 
oo  aa  JM.RDH 

I.J.DP8I(I.J).RDP8I(I,J) 

WRITE!*, 0818) 

jmiTE!*.*OIS)  <RV(  I).  1*1,  IAC) 

WRITE!  *.*0 17) 
oo  as  1*1, ROR 
JAC* JACC( t) 

WRITE!  S.SS 18)  I 

WRITE! S.SS 19)  !DV! I, J) ,J* l.JAC) 

CALL  COVAXB  FOR  CALCOLATIOR  OF  THE  BASIC  ODARTITIES  OF  THE  OOVARI- 
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3# 


CALL  COVAXSI -l.COV) 

CALCULATE  ALL  RECESS ARY  C0VARIARCE8  AT  ALL  GRID  POIRTB 


CALCULATE  COV(T.T)  AT  ALL  CRID  POIRTB 


I AC la  I 
JAC1*  1 
IACU-# 

IRK*# 

IRK-  IRK*  I 
IACL-IACU 
IACU- I ACC(  IRK) 

I ACL* IACL+1 
JAC* JACCI IRK) 
KI(6>*1 
KI( 7) * 1 

CALL  COVAXSI  # , COV) 
DO  4 I* I ACL, IACU 
CR( 2)  *RV(  I) -RE 
CRI3)-CRI2> 

DO  4 J-I.JAC 
CR(  I ) *DV<  IRK.  J) 
CALL  COVAXSI  l.COV) 
C< 1 . J. I. l)*COV 


CALCULATE  COV( DR(  T) . DR' ( T) ) AT  ALL  GRID  POIRTB 


KI(6)*2 

KI(7)*2 

CALL  COVAXSI  # , COV) 

DO  8 I* I ACL, IACU 
CR(  3)  *RV(  I) -RE 
CR(  3)  *CR(  2) 

DO  8 J* I , JAC 
CR(  1 ) * DV(  IRK, J) 

call  covaxsi i,cov> 

Cl 1.J.2,  l)*COV/RV(  D^RVI  I) 


CALCULATE  COV(  DTI  T) , DT( T) ) AT  ALL  GRID  POIRTB 


KII61-12 
KK7)  • 1 

CALL  COVAXSI • , COV) 
DO  6 I-IACL, IACU 
CR<  2)  *RV(  D-RE 
CR(  3)  *CR<  2) 

DO  6 J-l.JAC 
CR<  1)*DV( IRK, J) 
CALL  COVAXSI l.COV) 
C<  I, J.3.  D-COV 


CALCULATE  COV( DELTA  G. DELTA  G)  AT  ALL  GRID  POIRTB 


KII6)*3 

Kl(7)*3 

CALL  COVAXSI  • , COV) 

DO  7 I-IACL. IACU 
CR(  2)  *RV(  I)-RE 
CR(  3)  -CR(2) 

DO  7 J-l.JAC 
CR(  I ) »DV(  IRK,  J) 

CALL  COVAXSI l.COV) 

Cl  I.J.4.  D-COV/RVI  D/’RVI  I) 


CALCULATE  FIRST  ORDER  DERIVATIVES  AT  ALL  GRID  POIRTB 


KII6)-6 
KII7)-  1 

CALL  COVAXSI#, COV) 
DO  6 J- I , JAC, JAC I 


- ^ ^ _ .xt,  i i 


B 


* 


IB 


II 


ia 


13 


14 


CM  1>>W<  IRK,  J) 

52  • l-IACL. I AGO 

CR<  2)  *RV(  D-RE 
CRI  3)  «CRI  3) 

CALL  COVAWI  I.COV) 

C< I.J, l.3)-COV 

KII4MIB 

KKT)>a 

CALL  COVAXBIB.COV) 

00  * J-I.JAC.JAC1 
CR( l)»DV< IRK, J) 

DO  4 !■ I ACL. I ACU 
CRI2)«RV<  I) -RE 
CRI3)*CRI2) 

CALL  COVAXBI 1 ,COV) 

C(  I, J,a.3)>C0V/RV( D/HVI I) 

kk4>*  ia 

KII7)>4 

CALL  COVAXBIB.COV) 

DO  IB  J-t.JAC.JACl 
CRI 1)«DVI IRK,  J) 

DO  IB  I* IACL, I ACU 
CR(  3) *RV(  I) -RE 
CRI3)>CRI2) 

CALL  C0VAX8( I.COV) 

Cl I.J.3,8)*C0V 

KI(«)>B 

KII7>-3 

CALL  COVAIBIB.COV) 

DO  II  J-l.JAC.JACI 
CR< l) *DV< IRK, J) 

DO  II  I- IACL, I ACU 
CRI2)«RV<  I) -RE 
CRI3)>CRta> 

CALL  COVAXBI  I.COV) 

C<  I . J , 4, 3)  ■COV/'RVI  D/RVI  I) 
Kii«)>a 

KI<7)«  I 

CALL  COVAXBI  B, COV) 

DO  12  !■ IACL, IACO, I AC I 
CRI2)*RVI  D-RE 
CR(  3)  *CR(  2) 

S«RB2/RV(  D/'RVI  I) 

DO  12  J* 1, JAC 
CRI  1)»DV(  IRK, J) 

CALL  COVAJBI  I.COV) 

C< I.J. 1.2>*COV/8 

KII4M3 

KK7XB 

CALL  COVAXBI  4 , COV) 

DO  13  !■ IACL, I ACU, 1AC1 
CRI2)*RVI  D-RE 
CRI3)«CRI2> 

DO  13  J«I.JAC 
CR< I ) »DVt IRK, J) 

CALL  COVAXBI  I.COV) 
C(I,J,2,a)-COV/RB2 
Kit  6) *12 
Ki(7)*a 

CALL  COVAXBI  B. COV) 

DO  14  I* IACL, IACU, I AC I 
CR(2)>RV(  D-RE 
CRI3)*CRta) 

8-RB2/RVI  D/RVI  I) 

DO  14  J* 1 , JAC 
CRII)*DVI1RK,J) 

CALL  COVAXBI I.COV) 

Cl  I, J,3,a>*C0V/8 

KII4M4 

KI(7)*3 

CALL  COVAXBI  B , COV) 

DO  18  I* IACL, IACO, I AC I 
CRI2)*RVI  D-RE 


j — I MM  ■ 
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CR(3)=CR<2) 

8«RB2/RV(  I)/RV<  I) 

DO  18  J= 1 , JAC 
CR(  1 ) »DV(  IRK,  J) 

CALL  COVAXSd.COV) 

18  Cd,J,4,2)»COV/RB2 

C 

C CALCULATE  SECORD  ORDER  HIRED  DERIVATIVES  AT  ALL  CRID  POIRTB 

C 

DO  16  I* IACL, IACU, IAC! 

CR(2)*RV(  I) -RE 
CR(  3)  »CR<  2) 

S=RB2/RV<  I)/RV<  I) 

DO  16  J* 1 , JAC, JAC 1 
CR(  1)  = DV( IRK, J) 

KI<  6) *2 
KI<7)«6 

CALL  COVAXS(O.COV) 

CALL  COVAXSl l.COV) 

C<  I,J, 1.4)=C0V/S 

KI<6)>8 

KI(?>  * 10 

CALL  COVAXS(O.COV) 

CALL  COVAXSl  1 , COV) 

C< I,J,2,4)*C0V/RB2 

KI(6>=16 

KK7)  * 12 

CALL  COVAXS(O.COV) 

CALL  COVAXSC l.COV) 

Cd,J,3,4)=COV/8 
Kl( 6) *4 
KI( 7) *8 

CALL  COVAXS(O.COV) 

CALL  COVAXSt l.COV) 

C<  I , J , 4, 4) =C0V/RB2 

16  CONTINUE 

IF(  IRK.LT.NDR)  GOTO  36 
DO  17  1=1. IACU 

17  RV( I)=RB2/RV( I)/RV( I) 

C 

C CALCULATE  POLYNOMIAL  COEFFICIENTS  FOR  ALL  ELEMENTS 

C 

DO  26  11=1, NDR 
NDN=NDPS( II) 

NDH=JACCI( II.NDN) 

DO  26  J= 1 , NDH 

DY* DV(  II, J+1)-DV(  II, J) 

NNN=IACC( II)-1 

IFdl.CT.  1)  NNN=NNN-IACC<  1 1—  1 ) 

DO  26  NN= 1 , NNN 
I = NN 

IFdl.CT.  1)  1=1+  IACC(II-l) 

DX=RV(  I+1)-RV(  I) 

26  CALL  BILDEC 
CALL  TEST 
STOP 
END 


I 


BLOCK  DATA 

IMPLICIT  REAL*8(A-H.0-Z) 

COMMON  /C0MDAT/D1 , D32,  D2, D3 , Dl# 

COMMON  /C0H6/K7< 14) , K9( 14) , K1 1(14) , KIS( 14) , K19< 14) , K21( 14) , 
*K23( 14),  Cll( 14)  ’ 

DATA  D1,D32,D2,D3,D16  /I .000, 1 .8D6.2.6D6.3.6D6, 1 .60-14/ 

DATA  K7/8*6, 6* 1,3*2/,  K19/1 ,4*0, 2*1 ,7*0/,  K21/6, 1,3*6, 6*1, 3*2/,  K2 
*3/6*6, 1,6, 1,6, 1,6, 1,2/,  Cl  1/1 .6D6.-1.6D9, 1.608,2*1. 6D9,  2*266264.8 
*6606, 2*-1.6D9, 2*1. 609.3*1.609/,  K»/S*1,2,3,2,S,2.3,2,2,3/,  Kll/U* 
*6, 2, 3, 3/.  K13/1 1*1 ,2,3,3/ 

END 
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subroutine  bildbc 

0CCCCCCCCCCCCCCCCCC0C0CCCCCCCCCCCCCCGGC0CGG0C<XXXXX2C0CCCCCCCCCCCCCCCCCCC 

c c 

C GENERATION  or  THE  BICUBIC  POLYNOMIAL  COEFFICIENT  VECTORS  TOR  THE  C 

C ELEMENT  WITH  LOVER  LEFT  COBBER  INDICES  M.R  . C 

C C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCGCCCCCCGOCCCCCCCCCCCCCCCCCCCCCCCCCC 


********** 


E*  1 .4 


N. 

BX. 


IK 


BY 


********** 


.E>,  K»  1 • 


INPUT 

4-D1H. 
<K»  I) , 
( K«3)  , 


14 


********** 


TRARSrER  IR  COMM  « D.  DX.  BY.  N.  R,  Dl,  DBS.  03.  09 


C 
C 
C 

ARRAY  CONSIST IRC  OF  FURCTIOR  VALUESC 
I.  DER.  IR  X (K»3>.  I.  DER.  IR  Y C 
3.  OER.  IR  XY  (K»4>  C 

INDICES  OF  THE  ELEMENT'S  LOVER  LETT  CORN EMC 
GRID  DISTANCES  OF  TVS  CORRESPONDING  RLE-  C 
RENT  C 

C 

OUTPUT  **********  c 

C 
c 
c 

C 

c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccaxxxxxxxxxxxxxccc 

IIVLICIT  REAL*8( A-H.O-Z) 

COmnR  /GRID/D<94.94.4.  I4>/COMa/DUHMYia>,  DX.  DY.  M.  N ^COHDAT^Dl , 

« ooa.  oa.  ds 

DIMENSION  CC( 4.4) .C( 14) ,A(  14) 

EQUIVALENCE  ICCU  . I)  ,C(  I)  ) 

DO  44  J* I ,4 
DO  84  KM  ,3.2 
KfLOiN 
El  l"K+l 
DO  84  L* I .3.2 
Ll-Lsa+R 
LI l*L4 | 

CCtK.L)>DtKl.Ll.J.l> 

CCt El  I ,L) *Dt El ,Ll , J ,a>*DX 
CC(  K.Lll)*D<KI,LI , J,3)*DY 
CC<  El  I , LI  I >aDt  El ,LI . J,4)*DX*DY 
C1«D2*(C(»)-C<  I ) > 

ca>  D32*C 1 

C3*  D3*( C( I4)-C(2) ) 

C4»  D32*C3 
C8-Da*(C< 1 1 ) — C< 3) > 

C4>ma»c8 
C7«D2*<C< ia>-C(4)> 
ca» D32*C7 
C4»C( I3)4C<  8> 

CI4»C4+C< 8) 

CII>C( 14>+Cf4> 

Cia>C114C(4) 

C13«C<  l8)*CtT) 

CI4>Cl9+CtT) 

CI8>C(  I4)4C(B> 

CI4"CtS+C(S> 

A(  D-CII) 

Ata)«c<l> 

A(8)aC(8) 

At4)aCt4) 

At4)>C2-CI4 
A<  I3)>C9-C1 
A(  l4)aC4-ClB 
At  I4)*CI  1-C8 

At4>aCt4)4Cta)-M*tCt8)-Ct  1>) 

A<3)>C<3)-C<a>-C(  I ) -A(  4) 

AtB>aCtS>*Ct4>-Da*tCtY)-Ct8)) 

At7>aC(7)-Ct4)-CtB)-AtS) 

At 121 *00-014441 14)-D0*(C4-C14-AI4>) 

At  ID aG4-CI4-At4)-At  I4)-At  IS) 

At l4)*Da»tC8-Cl84At I9))-C74C18+A< 14) 

At IB)aCI3-C8-At 18) -At I4)-At 14) 

DO  41  11*8. 14 


~ * 


• li&taMfti  ...  .. 
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41  D(M.N.J.  II)»A(  II) 
4*  CONTINUE 
RETURN 
END 


SUBROUTINE  TEST 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccctxccccc 

c 

THIS  PROGRAM  IS  DESIGNED  AS  A USER  EXAMPLE  FOR  COVAPP  AND  PRO- 
VIDES, MOREOVER.  A COMPARISON  OP  COVAPP  AND  COVAXN,  A SLIGHTLY 
MODIFIED  VERSION  OF  THE  ORICINAL  SUBROUTINE  CO VAX  (TSCHERNINC, 

1474) , 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


THE  NAMES  OF  PARAMETERS  AND  VARIABLES  BETWEEN  THE  COMMENT  CARDS 
Cl  1 11  - CIIII  ARE  IDENTICAL  WITH  THOSE  USED  IN  COVAPP,  BETWEEN 
C2222  - C2222  WITH  THOSE  USED  IN  CO V NET  AND  COVAXN 

ADDITIONAL  ARRAYS  USED  FOR  COMPARISON  PURPOSES  : 

ICH( . ) ...  VECTOR  USED  FOR  REORDERING  THE  COVAXN  - 

COVARIANCES  (8.4. 14, ID  > (14.11,8,4) 

SEC< . , . ) . RES( . , . ) ...  ARRAYS  USED  TO  STORE  ALL  COVARIANCES  CAL- 
CULATED BY  COVAPP  AND  COVAXN 

CT( . ) ...  VECTOR  USED  FOR  COORDINATE  TRANSFER  TO 

COVAPP 


IF  COVNET  AND  TEST  ARE 
•C  GET  COVARIANCE 

HAS  TO  BE  REPLACED  BY 
CRID  PARAMETERS  SUCH. 
LABELED  COMMON  BLOCKS 

********** 

XLATP.  XLONP,  HP  ... 

XLATO,  XLONQ,  HO  . . . 


********** 


RUN  SEPARATELY.  THE  COMMENT  STATEMENT 
AND  GRID  PARAMETERS' 

A READ  STATEMENT  FOR  THE  COVARIANCE  AND 
THAT  ALL  QUANTITIES  ENCOUNTED  IN  THE 
/GRID/  AND  /CR1 PAR/  ARE  TRANSFERRED  TO  TE8TC 

C 

INPUT  **********  C 

C 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


CEODETIC  LATITUDE. 
OF  THE  POINT  P 
CEODETIC  LATITUDE, 
OF  THE  POINT  Q 


LONGITUDE  AND  HEIGHT 
LONGITUDE  AND  HEIGHT 


OUTPUT 


********** 


ALL  TYPES  OF  COVARIANCES  (1-14.1-14)  CALCULATED  WITH  COVAPP  AND 
WITH  COVAXN  AS  WELL  AS  ITS  ABSOLUTE  DIFFERENCES. 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
IMPLICIT  REAL*8<  A-H.O-Z) 

C 

Cl  Ml  ARRAYS  AND  DATA  NECESSARY  FOR  COVARIANCE  CALCULATIONS  USING  COVAPP 

C 

COMMON  /CRID/C(34,54.4, 14)  /CRIPAR/RB2,  S4<3>,  DS(3>,  Sl<3>.  DP8I( 

* 3,3).  PSIS(3,4),  PSIM3).  PS<3,84).  RE.  1ACC(3),  J ACC  1(3,8) 

* . NDR  /CONl/CTI  12) 

DATA  RHOC.GM  /8. 72437748 1D1 . 3. 48D14/ 

C 

CM  1 1 
C 

C2222  ARRAYS  NECESSARY  FOR  COVARIANCE  CALCULATIONS  USING  COVAXN 

C 

LOGICAL  LOCAL 

COMMON  /CNCOV/CK  12) , CR(8t).  SIGMA4(344>,  SICNA(344).  KM28).  Nl. 

* LOCAL 
C 

C2222 
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C3333  ARRAYS  ARD  DATA  USED  FOR  COMPARISON  PURPOSES 

C 

DIMENSION  I CHI  14) 

DIHEN8I0N  SEC( 14, 14) , RES( 14. 14) 

DATA  ICH/I.2.3,4,8.4.7. 14. 1 1 ,8.4. 12, 13. 14/ 

0 

C3333 

C 

3403  FORMAT!  1615) 

5414  FORMAT! 3D24. 14) 

6444  FORMAT!  1H1.4X.  'GEODETIC  COORDINATES  OF  P AND  0',//7X, ' I LATI 

*TUDE!P)  LONGITUDE!  P)  HEIGHT! P)  LATITUDE! ft)  LONGITUDE 

*!0)  HEIGHT!  O)  V/) 

6414  FORMAT!  1H  . 4X. 14. 2! 2F15.6 ,FI2. t ,5» > 

6411  ^FORMAT!  1H4.  18X.  • I* , 14X.  '2' . 14X.  ’3’ . 14X.  *4* . 14X.  ’5* . 14X,  '6' . 14X,  *7* 

6412  FORMAT!  1H  ,4X, I2.7F15.8) 

6413  FORMAT!  1H4, 4X, ’ MATRIX  OF  COVARIANCES  CALCULATED  WITH  COVAPP'/V) 

6414  FORMAT! 1H4, 18X. *B’ . 14X, , 13X, • 14* , 13X. • II* , 13X. * 12* , 13X, * 18* , 13X 
V, * 14’//) 


6415 

641S 

C 

C 

C 

C 

C 

C 


FORMAT! 1H4.4X, 'MATRIX  OF  COVARIANCES  CALCULATED  WITH  COVAXR'/V) 
FORMAT!  1H4. 14X. 'RELATIVE  ERRORS’ /V) 

GET  COVARIANCE  AND  GRID  PARAMETERS 


READ  AND  ECHO  THE  CE3DETIC  COORDINATES  OF  P AND  O 


NRDAT* 1 

444  READ! 5.5414. END-44)  XLATP , XLONP , HP , XLATO, XLONO, HQ 
WRI TE( 6 6004) 

WRITE! 6 6414)  NRDAT, XLATP, XLONP, HP, XLATO, XLONO, HD 
XLATP*  XLATP/ RHOC 
XLONP*  XLONP/RBOC 
XLATO*  XLATO/ RHOG 
XLONO*  XLONO/ RHOC 
C 

C CALCULATE  CR!  1) CR!  12) 

C 


CR!  2)  *HP 
CR!  3) *HO 

CR!4)*D6 IN! XLATP) 

CR!  5) *DSIN! XLATO) 

CR!  6) -DCOS! XLATP) 

CR!  7) * DC06! XLATO) 

CR! 8) * DS I N! XLONO- XLONP) 

CR!  9 ) * DC OS! XLONO- XLONP) 
P8I*CR!4)*CR!  5) +CR!  6)  *CR!  7)  *CR!  9) 
CR!  1)*PSI 

CR!  12)  *DARCOS!  PSD 

CR!  14>*CN/!CR!2)*RE)/!CR!2)+RE) 

CR!  11)  *CJH  CR!  3)  ♦RE)  /!  CR! 3)+RE) 


DO  744  1*1,  12 
744  CT! I) *CR( I) 

C 

C CALCULATE  ALL  DIFFERENT  TYPES  OF  COVARIANCES  US INC  COVAPP 

C 

DO  144  J*l, 14 
DO  144  K*  1 . 14 
CALL  COVAPP! J.K.COV) 

RES! J.K)-COV 
SEC!  J , K)  * RES!  J , K) 

144  CONTINUE 
C 

C PRINT  ALL  COVARIANCES 

C 

WRITE!  6. 64 13) 

WRITE! 6,6411) 

DO  344  J*  1 , 14 

344  WRITE! 6, 64 12)  J , ! RES!  J, K)  . K*  1 ,7) 

WRITE! 6,6414) 

DO  444  J*l,  14 

444  WRITE! 6, 64 12)  J,!RE8! J.K) ,K*8, 14) 
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C 

C CALCULATE  ALL  DIFFERENT  TYPES  OF  COVARIANCES  US INC  COVAXN 
C 

CALL  COVAXN! -1, COV) 

DO  600  1=1,  14 
KI( 6)  = I 
DO  600  J= 1 , 14 
KI(7)=J 

CALL  COVAXN!  O, COV) 

CALL  COVAXN ( l.COV) 

C 

C REORDERING  OF  COVARIANCES  (8,9,10,11)  > (10,11,8,9) 

C 

II=ICH( I) 

JJ= ICH( J) 

RES( II,JJ)=COV 
600  CONTINUE 

WR1TE(6.6015) 

WRITE(6,6011) 

DO  602  J= 1 , 14 

602  WRITE! 6, 60 12)  J,(RES( J.K) ,K* 1,7) 

WRITE!  6, 60 14) 

DO  603  J= 1, 14 

603  WRITE(6,6012)  J,  (RES( J.K) ,K=8, 14) 

DO  567  1=1,  14 

DO  567  J=l,  14 

SEC( I , J) = 1 . ODO-DABS( SEC( I, J)/RES( I, J)) 

567  SEC( I , J) =DABS(  SEC(  I , J) ) 

WRITE(6,6018> 

WRITE!  6,6011) 

DO  500  J= 1 , 14 

500  WRITE! 6,6012)  J, (SEC! J.K) ,K= 1 ,7) 

WRITE! 6, 60 14) 

DO  501  J=  1,  14 

501  WRITE! 6 , 6012)  J,  ( SEC!  J, K) , K*8, 14) 

NRDAT=  NRDAT+ 1 

GOTO  444 

40  RETURN  ~i' 

END 


SUBROUTINE  COVAPP! KR1 , KR2 , COV) 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c c 

C C O V A P P COMPUTES  APPROXIMATIONS  OF  A COVARIANCE  BETWEEN  TWO  C 
C QUANTITIES  RELATED  TO  THE  DISTURBING  POTENTIAL  OF  C 

C THE  EARTH.  C 

C C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C c 

C IT  IS  ESSENTIAL  TO  STATE  TP  T THIS  PROCRAH  ALONE  18  NO  ALTERNA-  C 

C TIVE  TO  THE  PARENT  - SUBRC  . INE  C O V A X ( C.C.  TDCHERN INC,  C 

C 1976) . COVAPP  IS  A DIFFERENTIATION  - INTERPOLATION  SUBROUTINE  C 

C WHICH  IS  ABLE  TO  PROVIDE  THE  SATE  KIND  OP  O0VARIANCBS  AS  COV  AX  C 

C DOES  WITH  ARBITRARY  HICH  ACCURACY  FOR  A MUCH  LOWER  PN'CE  IN  C 

C TERMS  OF  CPU  - TIME.  C 

C IT  IS  DESIGNED  FOR  LARGE  SCALE  APPLICATIONS  OP  COLLOCATION  IN-  C 

C VOLVING  THE  CALCULATION  OF  A LARGE  NUMER  OP  COVARIANCES.  C 

C THE  BASIC  NETWORK  OF  COVARIANCES  IS  GENERATED  BYCOVNET  C 

C USING  COVAXS,  A SIMPLIFIED  VERSION  OF  COV AX.  THIS  NET  HAS  TO  BE  C 

C GENERATED  ONCE  FOR  ALL.  AS  SOON  AS  IT  EXISTS,  COVAPP  CAN  BE  ACTI-C 

C VATED.  C 

C C 

C FOR  REFERENCE  SEE  (H.  SUENKEL,  1978).  C 

C C 


oftfioooooonfiooofiftoooorinnoonnfto  nflnftoonnnnnonnflnoooonnnonnfion 


Al( 

A2( 

A3< 

A4< 


********** 


POLYNOMIAL  COEFFICIENTS  OF  COV(T.T) 

• COVCDR(T)  , DR’  < T)  ) 

• COV(  DT<  T) , DT<  T) ) 

• COV< DELTA  C, DELTA  C) 


INPUT 


********** 


KR1 , KR2 


IN  ORDER  TO  HAKE  COVAPP  AS  COMFORTABLE  AS  POSSIBLE,  WE  HAVE  HINI-C 
HIZED  ITS  DISTANCE  TO  COVAX  RESULTING  IN  PARTLY  IDENTICAL  STATE-  C 
RENTS.  C 

C 
C 
C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


CR<  1) 

cr<  a> , 

CRM)  , 
CRM)  , 

cr<  a> . 

CR(  IS)  , 
CR(  12) 


CR(3)  . 
CR<  8)  . 
CR(7)  . 
CR(9)  . 
CR(  11) 


INTEGERS  CODING  THE  KIND  OF  OPERATOR  TO 
BE  APPLIED  ON  THE  COVARIANCE  FUNCTION  AT 
THE  POINTS  P AND  0 

COSINE  OF  THE  SPHERICAL  DISTANCE  BETWEEN 
P AND  O 

HEIGHTS  OF  P AND  Q ( METERS) 

SINE  OF  THE  LATITUDE  OF  P AND  Q 
COSINE  OF  THE  LATITUDE  OF  P AND  Q 
SINE,  COSINE  OF  THE  LONGITUDE  DIFFERENCE 
REFERENCE  GRAVITY  AT  P AND  O 
SPHERICAL  DISTANCE  BETWEEN  P AND  ft 


THE  VECTOR  CR  IS  TRANSFERRED  IN  THE  LABELED  COMMON  BLOCK  /COH1/. 
ALL  ELEMENTS  CONTAINED  IN  THE  LABELED  COMMON  BLOCKS  /CRIPAR/  AND 
/GRID/  MUST  HAVE  BEEN  MADE  AVAILABLE  BEFORE  COVAPP  IS  CALLED. 

TO  ALL  ELEMENTS  CONTAINED  IN  THE  LABELED  COMMON  BLOCKS  /COMB/ 

AND  /COMDAT/  A VALUE  HAS  BEEN  ASSIGNED  IN  THE  BLOCKDATA  SUB- 
ROUTINE. THE  MEANING  OF  K7( . ) K23( . ) AND  Cl  K.)  18  I DENT  I -C 

CAL  TO  THAT  IN  COVAXN.  C 


RB2 


RE  . 
NDR 


S«( I) , I* 1 , NDR 
DS( I) , I-l.NDR 
SKI). 1*1, NDR 

DPSK1.J)  ... 
PSIS(I.J)  ... 
PSIKI)  ... 
PS(I.K)  ... 


SQUARE  OF  THE  RATIO  BETWEEN  THE  RADIUS  OF 
THE  BJERHAHMAR  SPHERE  AND  THE  MEAN  RADIUS 
OF  THE  EARTH  (RE) 

MEAN  RADIUS  OF  THE  EARTH 

NUMBER  OF  RANGES  IN  RADIAL  DIRECTION  ( ITS 
MAXIMUM  IS  3 ACCORDING  TO  THE  DIMENSION 
STATEMENTS) 

LOWER  LIMITS  OF  RADIAL  RANGES  IN  S 
GRID  DISTANCES  IN  8- DIRECT I ON 
UPPER  LIMITS  OF  RADIAL  RANGES  IN  S 

SPHERICAL  DISTANCE  CRID  8PACINGS 
SPHERICAL  DISTANCE  RANGE  STARTINC  VALUES 
SPHERICAL  DISTANCE  UPPER  LIMIT  VALUES 
EACH  ROW  VECTOR  CONTAINS  THE  COSINE  OF  THEC 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


SPHERICAL  DISTANCE  OF  THE  CRID  POINTS 


********** 


OUTPUT 


********** 


cov 


COMPUTED  COVARIANCE 


THE  PROGRAM  U8E8  THE  FUNCTION  B8FC, 
BLOCKDATA. 


THE  SUBROUTINES  GET  AND 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
IMPLICIT  REAL*8( A-H.O-Z) 

COMMON  /AA/A1( 16) , A2< 16) ,A3( 16) ,A4( 16) ,D(36) 

/CRIPAR/RB2,  86(3),  DS<3>,  81(3) , DPSIO.B),  PSI8<3,6), 

P81 1(3) , P8<  3,(16)  , RE.  IACC(3),  J ACC  I (3,0)  , NDR 
/C0N»/K7(  14)  , K»(  14)  , K1 1(  14)  , K13(  14)  , K19(  14)  , K21(  14)  , 
K23( 14) , Cl  1(14)  /COMl/CR(ia>  /COM2/X,  Y.  DX,  DY,  IS. 

JT  /CO® AT/D  1 , D32,  D2.  D3 

EQUIVALENCE  <T,CR< 1) ) , (HP,CR(2> ) , <HQ,CR(3) ) , (SP,CR(4) ) , (SQ.CRI 8) ) 
, (CP,CR<  6) ) . ( CQ, CR( 7) ) , ( SD,CR(  8) ) , ( CD, CR( 9) ) , ( PSI , CR< 12)), 

< SSS,  D(  ID)  , (88C,D(29)  > • (8C8,D(  14)  ) . (SCC, D(  32)  ) , ( CS8, D( 28)  ) 
. (CSC,D( 12) > , ( CCS, D( 16) ) , (CCC,D( 18) ) 

6664  FORMAT!  IH  , 4X,  * 80RRY,  POINT  LOCATION  OUTSIDE  LEGITIMATE  LIMITS*/) 
R*HP+RE 
RS-HQ+RE 
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KCl'KRI 

KG3-KR2 

C 

C CALCULATION  OF  THE  PROPER  UNITS  OP  THE  00  VARIANCES  <CR1S> 

C 

CR13-CI  t(KRl)*CU!KRa> 

IP!  K19!  KR1  > . EA.  1 ) CRia«CRia/CR( 10) 

IF(K19(KRa> .EA. 1)  CRia*CRia/CR< 11) 

IF(KaKKRI)  . EQ.O)  COTO  98 

cRia-cRia/R 

i r< kri.lt. ia>  goto  96 

CRia-CRia/R 

IF(KR1.E0. 14)  CRia«CRia/CR(6) 

94  IF<  K3S( KR1) .CE. I)  CRia*CRia/CR!6> 

98  IF!K31!KRa> .EA.O)  GOTO  94 

CRia«CRia/R8 
IF< KR3.LT. ia>  GOTO  97 
CRia*CRia/R8 

IF!  KR2.EA. 14)  CRia«CRia/CR<7) 

97  IP(K33(KRa) .CE. 1)  CR13>CR13^CR( 7) 

94  CONTINUE 

C 

C DETERMINATION  WHETHER  THE  ROLE  OF  R AND  RS  HAVE  TO  RE  INTERCHANGED 

C 

IFI KR2.CE. KR1 ) GOTO  I 

R88*RS 

RS*R 

R*R8S 

KCI-KR3 

KC3-KR1 

1 CONTINUE 
C 

C DETERMINATION  OF  THE  NUMBER  OF  HORIZONTAL  DERIVATIVES  Nl 

C 

NH*K7<  KR1 ) +K7( KR3) 

C 

C CALCULATION  OF  THE  ELEMENT  INDICES  AND  THE  POINT  POSITION  RELA- 
C TIVE  TO  IT 

C 

S*RB2/R/R8 

IF(8.LE.S9(  1)  . AND.S.GE.SKNDR)  > GOTO  3 

2 WRITE!  6. 6404) 

C0V>9.99D44 

RETURN 

3 1*9 

4 1*1+1 

IF(S. LT.S1(  I) > GOTO  8 
COTO  6 

8 IF(S.GT.80(  I+l>)  GOTO  3 

GOTO  4 

6 DX*D8( I) 

X*(S-S9< I) )/D8( I) 

IS*  DINT!  X) 

X-X-IS 

18*18+1 

IF!  I.GT.  1)  18*  18+ 1 ACC!  I- 1) 

IF!  P8I . LE. P81 1! I))  GOTO  7 
coto  a 

7 J*0 

8 J»J+1 

IF!  P8I . LE.  P8I8!  I.J+D)  GOTO  9 
coto  a 

9 JT*DINT! ! P8I-P8I8! I , J) )/DP8Ii I , J) )♦ 1 
IF! J.GT. 1)  JT-JT+JACCII I.J-l) 

Y-T-PS! I.JT) 

DY*P81  I.  JT+D-PS!  I.JT) 

Y-Y/DY 

C 

C CALCULATION  OF  THE  PARTIAL  DERIVATIVES  OF  000! PSI)  WITH  RESPECT  TO 
C THE  SPHERICAL  COORDINATES  OF  P AND  O (Cl 04),  IF  NECESSARY 

G 

C1*D1 


.j.  ... : - *»-■  «■» — 


38 


983 


!*£■■>  IM I,  Mil, m 

i-n<nu> 

J-KlHKRl) 

K-»(KR2) 

»K11(KR2) 

J1*K13(  K1U) 
ki-ki3<kr2) 

D<  1)*01 
C8»CP*80 
SC>SP*CQ 
8CC>SC*C0 

cc-cp*cot 

CC8«CC*SD 
C8C»C8*CD 
D( 2) "CS-8CC 
D<  3 > ■ CCS 
D(  7)  * SC- CSC 
D( 18) « -CCS 

irOH.CT.  1)  GOTO  980 
C1»D(  I+») 

GOTO  lOOi 
D(4)»-T 
8C8*8C*8D 
88*SP*8Q 
CCC-T-88 
SSC>SS*CD 
C8S*C8*8D 
0(5)>-SCS 
0(6)>-CCC 
D<  8)  “CC+SSC 
D(9)«-CSS 
D(  19) »-T 
D(3I)«-CCC 
HN>K6+6*M 
IJ«I+J 
C1>D( IJ+KJD 
IF(KH.CT.2)  GOTO  981 

WTO(  }^(J,,*D(B5+1)*ik«*«i-8) 

88*80 

D(  19)»-D(7) 

D(  17) —800 
D(  18)* CCS 
0(20) >-0(3) 

D(Sl)*~CC8 
D( 26) >-888 
D(27)«-C8C 
0(33) a-CG 8 
»«6*Hl-6 

ir(lfH.CT.3)  GOTO  982 

IjK  '♦*>*«•>■»«> 

0(22)  >T 

0(23) >SCS 

D(24>*CCC 

D(  28) »-CSS 

0<  39) «-C88 

D(34)«CCC 

DOS) -8C8 

D( 36) «CCC 

K»1»B+1 

N61-K+1 

UI>tt>I 

BJ-B6+J 

nti-m+i 

MU-HS+J 

KKl-Dt+1 

•SnaaS®  s>»e ***  '•>*»>-■>-»'.♦»<  >j> 

AmB\  I)*D(K61) 

B*D(  J)*D()KI) 


w I' 

981  888 


'■9X— 


— i 


C3-D(  I J ) *D(  K6 1 ) *D<  Ml)  +B*D(  K6 1 ) +D(  J)  *(  D(  N6 1 ) *D(  K6  1 ) +D(  Dtl ) *D(  I))+A 
*#D<N6J)+D(K6J)*D<  I)*D(N61)  1 * 

C4*A*B 
GOTO  1094 
CONTINUE 

IF(NH.CT.O)  KC2* < KG2-4) /2 
GOTO( 94 1.942. 943, 944, 945) .RC1 
CALL  CET(  1) 

GOTO< 11,12. 13, 14. 15) ,KC2 

COVARIANCES  INVOLVING  LE. 1 HORIZONTAL  DERIVATIVES 

COV-BSFC(AI,4.NH)*Cl 
COTO  1444 

C0V*-BSFC(A1, 1,NH)*8/R8*C1 
COTO  1444 

COV- ( BSFC<  A1 . 1 , NH) *0-B8FC( A1 . 4. NH) *D2> /RS*C1 
COTO  1444 
CALL  CET<  2) 

COV* - ( < BSFC( A 1 , 4 , NH) *D2+B8FC< A1 , 1 , NH) *S) /’R8-BSFC1 A2 , 4 , NH)  *R>  /R8*C 1 
GOTO  1444 
CALL  CET<  2) 

COV> ( BSFC( A2 . 4, NH) *R+B8FC( A 1 . 1 . NH) *8/R8) /'RSVC I 
GOTO  1444 
CALL  CET<  2) 

COTO(21 .22.23,24,25) , KG2 
CALL  CET<  I) 

C0V«-B8FC< A1 , 1 , 1)*8/R*C1 
COTO  1444 

COV*BSFC( A2.4.NH) *C1 
COTO  1440 
CALL  CET(  1) 

COV* ( -BSFCt A2 , 4 , NH) +BSFC< A1 , 1 , NH)  *D2*S/R/R8> *C 1 
COTO  1440 
CALL  GET!  1) 

COV*  ( BSFC(  A2 , 4, NH) *D2+<  BSFC( A1 , 1 , NH) *D2/R/RS-B8FC< A2. 1 , NH) > «8> /US 
**C1 

COTO  1444 

COV* -BSFC( A2 . 1 . NH) *S/RS*C 1 
GOTO  1444 

C0T0(31 ,32,33,34,35) , KGS 
CALL  CET< 1) 

COV*-(  BSFC(  A1 ,4, 1)*D2-BSFC(A1, 1. 1)*S)/R»C1 
COTO  1444 
CALL  CET(  1) 

CALL  GET( 2) 

C0V-(BSFC(A1, 1.1)*D2*S/R/RS-B8FC(A2.4. 1))*C1 
COTO  1444 
CALL  CET<  4) 

COV- B8FC( A4 , 4 . NH) *C 1 
COTO  1444 
CALL  GET(  4) 

COV*  BSFC( A4, 1 , NH)*8*C1/R8 
COTO  1440 
CALL  CET( 1) 

GALL  GET(  2) 

#COV* ( BSFC< A2 , 1 , NH) *8-< B8FC( A2 . 4 . NH) +BSFC( A 1 , 1 , NH) *S/R/RS) *02) /RS*C 
COTO  1444 

COT0141 ,42,43.44,45) , KG2 
CALL  CET< 1) 

CALL  CET( 2) 

COV-  <B8FC(A2.4, 1>*RS-<B8FC( Al. 1. 1)*8+B8FC< Al,4, 1)*DS)/'R)/R*C1 
GOTO  1444 
CALL  CET(  1) 

CALL  CETI2) 

C5S£“<B8^C,  *•  I>*8“<B8FC<A1,  1, 1>*8/R/R8+BSFC(A2,4, 1>)*D8)/'R*C1 

GOTO  1444 
CALL  CET(  4) 

COV*  BSFCIA4, 1, 1)*8/R*C1 
GOTO  1444 
C.Al-l-  CET(  4) 

COV* ( B8FC( A4.2.NH) *S+B8FC( A4, 1 , NH) )/R/RS*Cl*8 


GOTO  1M# 

48  CALL  CITI  I) 

CALL  GETO) 

cov<<B8rc<  A2,a,ra>*8-B8FC<Aa.  i,ni>-<BSFCtAa.4.iii)/soBrc<Ai.  i.rd/ 
*R/R8)  ) *8/R/RB*C  1 
GOTO  1444 

9M  call  err;  s> 

GOTO! 8 1.88,88, 88. 88> . 868 
81  CALL  CET(  1) 

COV-  (B8FC(  A2.4.  1 >KRS+B8FC< Al.  1. 1 ) *8SR> /IWC 1 
COTO  1444 
88  com  HUE 

COV— B8FC<A8. 1,1)*S/R4C1 
GOTO  1444 
88  CALL  CET< I) 

COV>  (B8FC(  A8, 1.  1)*8-<B8FC< Al,  I,  1)*8/R/RB+B6FC< AS.4,  l))tN)/1MI 
GOTO  1444 
88  COHTIRUE 

COV-  < B8PC(  A8 . 8 , 1IH)  *8*B8FC<  A8 , 1 , HE)  ) *8*C  1/R/R8 
GOTO  1444 
1498  C08TI HUE 
C 

C COV ARI AJICES  INVOLVING  LE.8  HORIZONTAL  DERIVATIVES 

C 

IFOGCl.GT.B)  GOTO  68 
CALL  CET< 1) 

CALL  GETO) 

GOTO! 118.818,318, 418. 8181.861 
118  COV-BSFC<A3,4,4>*Ca*B8FC<Al,4,l)*Cl 
GOTO  1444 

818  COV— <B8FC<  A3, 1,4)*C8+B8FC<A1, 1, 114C1148/R 

GOTO  1444 

318  C0V<<<B8FC<  A3, l.4)*9-B8FC<A3,4,4>*D8)*Ca*<B8FC<Al, 1. 1)*#-BBFC<A1,4 

*, I)*D3)*C1)/R 
COTO  1444 
418  CALL  CETO) 

C0V<<<B8FC<  A8, 4. 1 ) *1M>R8-B6FC<  Al  .4, 1>*D8-B8FC<  At , 1, 1)*8)*C1*<  BSFC<  A 
»3 . 8 . 41 *8*8-B8FC< A3 . 4 , 4) «D8) «C8) /R/R 
COTO  1444 
818  CALL  GET<  8) 

COV-  < < B8FC<  A8.4, 1) *R8+B8FC<  Al , 1 , 1)*8/R)«C1+<BBFC<  AS,8,4)*8+BSFC< A 
*3 , 1 , 4) *D8> *S/R*C8> /R 
GOTO  1444 

68  863- < 868-4) O 

861- <861-41 /8 
G0T0<64,61 ,68) ,861 
64  CALL  CET< 1) 

CALL  GETO) 

GOTO< 66,614,68) ,868 

66  C0V-B8FC<A3,4.4)*C8+B8FC(A1,4, 1)*C1 

COTO  1444 

68  COV-<  < B8FC<  A3, 1 ,4)*S-B8FC< A8,4,4)*D8)*CS+< BSFC< Al , 1 , 1)*8-B8FC< Al ,4 
*, 1)*D8)*C1)/R8 
GOTO  1444 

614  COV— <B8FC<A3,1,4)*C8+B8FC<A1, 1, t)*Cl)*8/RS 
GOTO  1444 
61  CALL  GET<8> 

GOTO! 88, 88. 814).  KGS 
88  CALL  GETO) 

Q0V-B8FC<A8,4, 1)«C1*(B8FC<A3,8,4)*8+B8FC< AS, 1,4) )*8/R/R8*C2 
COTO  1444 
814  CALL  CET< I) 

CALL  GETO) 

COV— <<B8FC<A3,8,4)*8-B8FC(A3, 1 ,4)  )«C8-B8FC<  Al , 1, 1 > *D8*C 1 ) *8^R/R8- 
*B8FC<  A2,4, I)*C1 
COTO  1444 
68  CALL  GET<  4) 

CALL  GETO) 

COV-  B8FC<  A4 , 4 , 1)*C14<BSFC<  A3.8,6)*8-BBFC<A3, 1 ,4)*D3*B8FC< A3,4,4)*D 
*8*D8>*8/R/RS»Ca 
COTO  1444 
1494  CALL  GET< 1) 

CALL  GETO) 


41 


C 

C 

c 

619 

1919 

819 

C 

C 

C 

1944 


1999 


COVARIANCES  INVOLVING  LE.3  HORIZONTAL  DERIVATIVES 


KC1-(HG1-4>S9 
coTo<6ia.aia. 1912) . not 

COV*B8FC( Al,6, 1>*C1+B8FC< A3,9,0>*Ca+B8FC<A3,9, t)*C3 
GOTO  1999 

COV*(<BSFC< Al, 1, 1>*CI*B8FC<A3, 1 ,9>*C9»B8FC< A3, 1, 1)*C3)*8-(BSFC<A1, 
*9,  l >*C1+B8FC<  A3,9,4)*C2-*-B8FC(  A3,9,  1>*C3>*D3>/R 
GOTO  1999 

C0V*-(B8FC(A1, 1, 1>*C1+B8FC(A3, 1 ,9)*C2+B8FC< A3, 1, 1>*C3)*8/R 
GOTO  1999 


COVARIANCES  INVOLVING  LE.4  HORIZONTAL  DERIVATIVES 


CALL  GET< 1) 

CALL  CETO) 
COV<BSFC<A1<  9,  1 

*4 

C0V>C0V*CR12 

RETURN 

END 


> *C l+BSFCI AS , 9 , 9> *C2+BSFC< A3 , 9. 1 ) *C3+BSFC< A8 , 9 , 2) *C 


SUBROUTINE  CETUO 

CCCCCCCCCCCCCCCCCCCCCCCCCOCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


C 

c 

C 

C 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


THIS  SUBROUTINE  MAPS  THE  POL YOU I AL  COEFFICIENTS  FOR  A SPECIFIC 
ELEMENT  < IS, JT)  ONTO  A VECTOR 


C 

C 

C 

C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


********** 


INPUT 


********** 


U<  . , . <K<  J> 


4-D1H.  ARRAY  OF  FUNCTION  VALUES  < J»l> , 
1.  DER.  IN  X (J>3),  1.  DER.  IN  V 

a.  DER.  IN  XT  <J-4>  FINN  EM 4 « 

HP  1 ...  COV<T,T), 

K*2  ...  COV< DR< T) , DR* < T) ) < 

R*3  ...  COV(  DTT  T)  , DT<  T)  ) < 

K*4  ...  COV< DELTA  G, DELTA  G> 


18,  JT 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


INDICES  OF  THE  ELEMENT’S  LOWER  LEFT  CORNERC 

C 

OUTPUT  **********  C 

C 

AH.),  A2< . ) , A3<  . ) , A4<  , ) ...  VECTOR  OF  POLYNOMIAL  COEFFICIENTS  C 
CORRESPONDING  TO  THE  ELEMENT  < 18,  JT) . C 


********** 


Alt  .) 
A2( . > 
A3<  . ) 
A4<  . ) 


POLYNOMIAL  COEFFICIENTS  OF  COV(T.T) 
COV<  DR(  T) , DR*  t T)  > 

COV<  DTX  T)  ,DT<T) ) 

COV< DELTA  G, DELTA  G) 


TRANSFER  IN  COMPS  > C.  IS,  JT,  Al,  A3,  AS.  A4 


C 

C 

C 

C 

C 

C 

C 


CCCCCCCCCCCCCCCCCCCCCCCCOOOQCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


1 

19 


IWLICIT  REAL*fi( A-B.O-Z) 

COMMON  /GR1D/C<99,S9,4, 14)  /COIB/DUNMYI 4) , IS,  JT 
» A9< 16) , A3< 16) , A4< 16) 

00T0< I. a. 3.4).  K 
DO  19  1*1,16 
Al<  I)«C( IS.JT.K,  I) 


/AA/A1(  16)  , 


DO  39  1*1,16 


42 


26 

A2<  I ) -C<  IS,  JT, K,  I) 
RETURN 

3 

DO  36  1-1,16 

38 

A3<  I)  -C<  IS.JT.K,  I) 
RETURN 

4 

DO  48  1-1,16 

43 

A4<  I)-C(  IS.JT.K,  I) 
RETURN 

END 

function  bsfcia, ix. m 

CCCCCCCCCCCCCCCCCCCCCCCCCOCCCCCCCCOCOOCCCCOCCCCCCCCCCCCCCCCCCGCCCCCCCCCC 

c 


THIN  FUNCTION  CALCULATES  ZERO.  FIRST  AND  SECOND  ORDER  DERIVA- 
DIVES  OP  A BICUBIC  POLYNOMIAL 


INPOT 


********** 


A(l«>  . 
IX,  IY 

x,  r .. 

DX.  DY 


VECTOR  OP  NORMALIZED  POLYNOMIAL  COEPP. 
ROWER  OP  X-  AND  Y-DERI  VAT  IVES 
NORMALIZED  COORDINATES  OP  IRE  CALCULATION 
POINT 

GRID  DISTANCES  IN  X-  AND  Y-DIRECTION 


********** 


OUTPUT 


********** 


I.  OR  2.  ORDER  DERIVATIVE 


THE  PROGRAM  USES  IRE  BLOCKDATA  SUBROUTINE. 
TRANSFER  IN  COHSMI  * X,  Y.DX.DY.Dl  .B32.D2.D8 


C 
C 
C 

a 

C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

IMPLICIT  REAL*8( A-H.O-Z) 

DIMENSION  A( 16) 

COMMON  /COHDAT/'Dl  ,DS2,D2,D3  /COie/X, Y, DX, DY 
0070(2,9) , IY 

BSFC>  A(  4)  ♦' Y*<  A<  8)  Y*(  A(  12)  +Y*A(  16)  ) ) 

GOTO(4.B> , IX 
C IX>6,  IY-6 

BSFC- A( 3) *Y*< A( 7) +Y*(  A(  1 1 ) + Y*A(  IS) ) > +B8FC*X 
BSFC-A<2)+Y*<A(6)*Y*<A(  16)+Y*A<  14>>)+B8FC*X 
BSFC-A( 1)+Y*(A(S)+Y*(A(»)+Y*A< 13) ) )+B8FC*X 
RETURN 

C IXil,  IY-6 

4 BSFC-A(S)+Y*<  A(7>+Y*(A( 11)+Y*A( 15) ) )+B8FC*X*D32 
B8FC- A( 2) +Y*(  A( 6) +Y*(  A( 16)+Y*A(  14) ) > +BSFC*X*D2 
BSFC-BSFC/PX 
RETURN 

C IXN2,  IY-6 

3 B8FC-A<3>*Y*(A<r>+Y*(A<  11)+Y*A<  18) ) >+B8FC*X*D3 

BSFC- B8FC*D2/DX/DX 
RETURN 

2 YD2»Y*D3 

YD82-Y*D32 

B8FC*A(8)+YD2*(A< 12)+YD02*A< 16)) 

GOTO<7,8) . IX 
C IX-4,  IY* 1 

B8FC>A<7)*YD2*<A<  1I)+YD82*A<  13)>+BSFC*X 
BSFC- A(  6) +YD2*(  A(  16M-YD32*A<  14))4B8fC*X 
BSFC-A(8)+YD2*<A<9)+YD32*A( 13))+B8FC*X 
B8FC*  BSFC/DY 


— Hi 


oononoflfiooooooofldoofio 


ix- 1,  nr- 1 

B8FC-A(7>+YD0*(  A<  II)+YD82*A<  IS) ) +MFC*]0*D82 
B8PC-AI6)  +YDS*<  A(  1*> +YD9a*A<  14) ) >MFC*K«Da 
B8FC- B8FC/DX/DY 
RETURH 
ix-a.  IY-1 

B8FC^A<7)4YD8*(A( ll)4YD82*A( 18) >+B0FG»X»DQ 
B8FC-  B8FOD2/DX/DX/DY 
RETURN 
IX-*.  IY-2 

B8FC«A(  13>+Xt<  A<  14>+)0ft<A<  I5)-»)0*A<  1«)>) 

B8FC- A< 9) +XM  A<  1*>+)0»<A<  IIH-XtAl  I8>  ) > ♦BSFC*Y»D8 

B8FC- B8FOKD2XDY/DY 

RETURN 


COVARIANCE  FUNCTION  PARAMETERS  < 


RBRE2 9 . 9996 179D+94 

A 0.4232844D+93  45 

KI( 3) 34 

KI(4) • 

KK3) 2 

N 12 


NUMBER  OF  RADIAL  RANGES  ( NDR)  ...  3 

I H9(I>  Hl(I)  DH(D 


1 


• .0 


200. • 


2 1M4M.4  1 39999.  • 7399. 9 

3 299999.9  299999. 9 19999.9 


NUMBER  OF  SPHERICAL  DISTANCE  RANGES 
I NDPS(I) 


1 

2 
3 


4 

3 

3 


1 


SPHERICAL  DISTANCE  GRID  VALUES  AND  NUMBER  OF  CORRESPONDING  INTERVALS 


1 

J 

DPSI( I,J> 

NDPSK 

1 

1 

0.999987 

19 

1 

2 

9.999291 

17 

1 

3 

4.999873 

13 

1 

4 

9.992999 

6 

2 

1 

9.992999 

6 

2 

2 

9.990727 

8 

2 

3 

9.917483 

19 

3 

1 

9.992999 

6 

3 

2 

9.998727 

19 

3 

3 

9.917483 

18 

GRID  POINT  VALUES  IN  RADIAL  DIRECTION 


9371999.9 

9372999.3 

9471999.9 

9991137.7 


9371299.9 

9372899.7 

947B399.9 

9911279.2 


9371499.9 

4373999.8 

9489929. 1 

9421491.8 


4371499. I 

4373291.9 

4493378.4 

4431494.4 


GRID  POINT  VALUES  IN  SPHERICAL  DISTANCE 


RADIAL  RANGE  NUMBER 


1 


1.4 

1. 

4.999997 

9.999988 

9.999948 

4.999883 


1.1 

l.l 

9.999994 

9.999989 
9.999939 
9.999799 


I. 

1.1 

9.999998 
4.999988 
4.999929 
9.999739 


1.4 

4.999999 

9.999994 

4.999983 

4.999918 

4.999448 


RADIAL  IIANCE  NUMBER  ...  2 


1.4 

9.999891 

9.992849 

4.998924 


9.999999 

4.999948 

4.994298 


4.999983 

4.998439 

4.987488 


RADIAL  RANGE  NUMBER 


3 


9.999391 

4.994822 

4.979294 


9. 

3. 

9.992844 

4.948924 


9. 

4.998489 

4.999248 

4.941242 


9. 

9.998138 

9.987488 

4.984398 


1 


4371899.1 

4373491.2 

4341 187.4 


4372999.2 

4373441.8 

4871999.9 


4372299.3 
4373891.7 

4881999.4 


4372499.4 

4374992.9 

4891948.8 


Ml 


1.4 

4.999999 
4.999993 
9.999978 
9.999997 
9.999888 


1.4 

9.999998 

9.999992 

9.999971 

9.999894 

9.999497 


1.4 

9. 

9 
9. 

9.999881 

4.999491 


1.4 
9.999997 

9.999989 
9.999987 
9.999867 


4.999942  4.999932 
4.998188  9.997844 
4.984898  9.981627 


4. 

9.996917 

9.978148 


9.999848 

4.994198 

9.974379 


4.999487 

9.994822 

4.979296 


9.999932 

4.997864 

9.984898 

9.981987 


4. 

4.994917 

9.981427 

9.948819 


4. 

4.996198 

9.978148 

9.939693 


9.999487 

9.998396 

4.974379 

4.933889 
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